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1. Introduction 


This report provides a detailed final summary of research conducted under 
NASA Grant NAG-1-657 entitled, “Reynolds Stress Closure in Jet Flows Using 
Wave Models.” 

The goal of the research has been to develop turbulence closure schemes, free 
of model constants, that predict the development of jet flows. In particular, the 
unsteady characteristics at the large scale were to be modeled and the models were 
to be applicable to jets of arbitrary geometry. In addition, the modeling scheme 
had to be computationally inexpensive. This eliminated the possibility of Direct 
Numerical Simulations or Large Eddy Simulations. 

The models developed under this program have achieved the goals of being free 
of model constants and being computationally inexpensive but the closure schemes 
were not applied to jet flows. Instead, the closure scheme was developed for a 
two-dimensional mixing layer. Extensions of the closure scheme to circular and 
non-circular jets flows is the subject of ongoing research activity under separate 
funding. However, the numerical procedures to apply the closure scheme to non- 
circular jets have been developed. This permits the local characteristics of the 
large scale structures to be calculated. In addition, the shock cell structure in 
non— circular jets, operating off— design, has been calculated. 

The modeling procedure is described in detail in the subsequent chapters. How- 
ever it may be summarized as the response to the following questions. 

1) If the time-averaged Bow properties are known, can the most iikeiy unsteady 
Bow Beld be deduced? 

2) If the time-averaged turbulent stresses associated with this unsteady Bow Beld 
are calculated, are they compatible with the original time-averaged velocity 
and temperature Belds? 

The answer to both these questions appears to be “yes” for the free shear flows 
examined. Such flows are dynamically unstable and are dominated by their most 
unstable, linear instabilities. The advantage of such an approach is that it models 
the unsteady flow field directly using a phenomenological model. It has no need to 



resort to model constants to close the higher order equations that are encountered 
in traditional models. However, the models have yet to be tested in more complex 
configurations. It is likely that models such as the one described in this report 
will act as guides for the modeling procedures in conventional Reynolds-averaged 
models. This will enable more physically realistic models to be incorporated in 
these closure schemes. 

The chapters of this report represent either publications or manuscripts sub- 
mitted for publication by the principal investigator and the research assistants. The 
three research assistants who have worked on this project each earned a Ph.D. de- 
gree. Many details of the research are given in their theses. The student names and 
their thesis titles are: 

Roy S. Baty, “Reynolds Stress Closure in Jet Flows Using Instability Wave 
Modeling,” Ph.D. thesis , Department of Aerospace Engineering, The Pennsyl- 
vania State University, 1989. 

Thonse R. S. Bhat, “Linear Models for the Shock Cell Structure of Super- 
sonic Jets with Noncircular Exit Geometries,” Ph.D. thesis, Department of 
Aerospace Engineering, The Pennsylvania State University, 1990. 

William W.— W. Liou, “Weakly Nonlinear Models for Turbulent Free Shear 
Flows,” Ph.D. thesis, Department of Aerospace Engineering, The Pennsylvania 
State University, 1990. 

The outline of this report is as follows. Chapter 3 contains a description of 
the weakly nonlinear turbulence model developed by Liou. An essential part of 
the application of such a closure scheme to general geometry jets is the solution of 
the local hydrodynamic stability equation for a given jet cross-section. Chapter 4 
describes the conformal mapping schemes used by Baty to map such geometries onto 
a simple computational domain. Chapter 5 describes Baty’s solution of the stability 
problem for circular, elliptic and rectangular geometries. In Chapter 6 Bhat’s use 
of linear models for the shock cell structure in non-circular jets is given. The 
Appendices contain reprints of papers also published during this study. Appendix 
I describes the instability of elliptic jets. Appendix II provides a technique for 



predicting the shock cell structure in non-circular jets using a vortex sheet model. 
Finally, Appendix III describes the resonant interaction between twin supersonic 

jets. 

Each of the component parts of this research program provide progress toward 
the prediction of the development of jets from arbitrary geometry nozzles including 
their unsteady turbulent structure at the large scale. It also provides the basis for 
the prediction of their near-field pressure fluctuations and their radiated noise. 
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WAVE MODELS FOR TURBULENT 
FREE SHEAR FLOWS 

W. W. Liou and P. J. Morris 
The Pennsylvania State University 
University Park, PA 16802 


ABSTRACT 

New predictive closure models for turbulent free shear flows are presented in this paper. 
Thev are based on an instability wave description of the dominant large-sca e s ruc ur 
in these flows using a quasi-linear theory. Three models have been developed to > stud* 
the structural dynamics of turbulent motions of different scales m free s ear _ , 

local characteristics of the large-scale motions are described using linear theory, 
amplitude is determined from an energy integral analysis. The mo e s ave ee 
to the study of an incompressible free mixing layer. In all cases, predictions are ma e or 
development of the mean flow field. In the last model , predictions of the 
motion of the large-scale structure of the mixing region are made. The predictions 
good agreement with experimental observations. 

I.INTRODU CTION 

Though the presence and importance of large-scale coherent structures to the mixing 
process in free shear flows has been recognized for many years turbulence mo d that 
incorporate these observations have been very limited. The use of direct numerical or large 
eddy simulations provide a detailed prediction of the large-scale motions in low h * 
Reynolds number turbulent flows respectively. But these predictions are computation } 
expensive and are still limited in general to simple boundary conditions. The prese 
model makes use of experimental observations in excited turbulent flows or conditional 
sampling measurements to provide a simple model of the large scale motions which is 
computationally inexpensive. 

Most current turbulent flow calculations for practical applications use the lon & tir ° € ' 
averaged Navier-Stokes equations. Turbulence models are needed to evaluate the unkno 
correlation terms, the Reynolds stresses, that appear when the statistical averaging process 
is applied to the nonlinear convective terms in the equations. This is the dosure Problem. 
There are closure models of various orders that have been proposed. These models are 
usually based on the notion that the high-order moments of fluctuations can be repre- 
sented reasonably well as functionals of moments of lower order. Work m this regar 
voluminous and will not be elaborated on here. Some models are quite success u an 
have become verv popular in engineering flow calculations. However, they usua y in ^° 
a large number of model constants determined by comparison with experimental da . 
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Thus, these models are not entirely predictive but, in some ways, represent a sophisticated 
correlation of experimental data. 

The present models are based on observations of large-scale coherent structures in free 
mixing layers. Brown and Roshko (1974), among others, observed that these orderly 
motions dominate the dynamics and the structure of free shear flows like wakes, jets and 
mixing layers. The structures appear in both low- and high-speed flows. They have also 
been observed in many flow geometries. 

This paper is concerned with new, predictive turbulence models for free shear flows. The 
models simulate the propagating large-scale structures as spatially travelling instability 
waves. In this paper, we focus on the validation of the wave models as well as a determina- 
tion of their limitations. Predictions are made for a two-dimensional incompressible free 
mixing layer. This will provide guidelines for applications of the models to more complex 
configurations. 


II. THE WAVE MODELS 

The wave models presented here are used to make a direct calculation of the large 
scale, characteristic structures. The fundamental idea is that the large-scale structures 
may be modeled using a quasi-linear theory. The local characteristics of these structures 
may be described by linear instability theory. This has been demonstrated by the ex- 
periments of Gaster, Kit and Wygnanski (1985) and Petersen and Samet (1988), among 
others. In their experiments they compared predictions of the amplitude and phase of 
the axial velocity fluctuations, based on linear stability theory, with phase-averaged mea- 
surements in an excited shear layer and a jet. The agreement between predictions and 
experiment was very good though only normalized distributions of amplitude and phase, 
not the absolute amplitude, were predicted. This close agreement between the predic- 
tions of linear stability theory and the properties of the large-scale coherent sturctures has 
formed the basis for theories of turbulent mixing and supersonic jet noise generation and 
radiation. For example, Tam and Morris (1980) and Tam and Burton (1984) predicted 
the noise radiation from instability waves in supersonic shear layers and jets and obtained 
very good agreement with experiment. Their analyses showed that the behavior of the 
large-scale disturbances could be modeled satisfactorily using a quasi-linear theory, even 
though the waves were not infinitesimal in magnitude. However, an important element of 
these calculations, the velocity profiles of the mean flow, that are needed for the linear 
stability calculations, are obtained from experiments. Their approaches provide a closure, 
but are not predictive. The models proposed here establish a complete closure model us- 
ing a simple quasi-linear theory for the large-scale motion. In the present model both 
the mean flow and the time-dependent turbulent motions at the large-scale are obtained 
simultaneously as a solution. 

Il.a Analysis 
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Turbulent flows generally contain broadband fluctuations and are traditionally treated 
as random processes. These processes are then described bv statistical averages^ T e 
observation and justification of the occurrence of orderly and coherent large-scale turbulen 
fluctuations in turbulent free flows encourage an interpretation of these turbulent flows 
using 21 quasi-deterministic description. 

In light of the existence of coherent structures, it is appropriate to decompose an 
instantaneous flow property into three parts. That is, 

u = f, + h + u (1) 


The fluctuation with respect to the mean quantity, F ,, is separated into two components, 
one representing the dominant large-scale motion, /*, and the other representing the sma - 
scale fluctuation, /'. The mean flow component is obtained from a long time-average ol 
its instantaneous value and is given by 


u = 


Fi = 



( 2 ) 


A short time-average may be defined by 


< fi > = Fi 



(3) 


where T 7 is much smaller than T u but much larger than the characteristic time scale o 
the background small-scale fluctuation. The mean flow represents an average over many 
realizations for a long period of time and thus is the profile that is most probably seen by the 
large-scale structures which occur randomly in space and time. In this approach, mean 
flow properties and large-scale fluctuations can be obtained explicitly. The small-scale 
turbulence which provides additional mixing at smaller scales compared to the mean an 
the large-scale motions is treated separately. The governing equations for the mean flow can 
be derived bv long time-averaging the Navier- Stokes equations. Equations governing t e 
large-scale fluctuations can then be obtained by subtracting the resulting mean equations 
from the short time-averaged Navier-Stokes equations. The governing equations for the 


mean flow' are 


dUj 

dxi 


= 0 


rr dUi d , tt\ _ dP 

v ’aT, + a^ {u<u ’ • r U( “' ) - “■ 


1 d 2 Ui 


(4 .b) 

(4. a) 


dxi Re dzjdxj 

where the interactions between motions of disparate scales are assumed to be negligible. 
The equations can be simplified further by applying the boundar 3 ’-laver approximatio 
For two-dimensional flow's, the resulting equations are 


8U dV_ 
dx ' dy 


0. 


(5.a) 
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(5.6) 


u™ + vg 

d x dy 




+ 


d_ 

dy 


(xi'v 1 ) 


1 d 2 U 
Re dy 2 


In order to close the governing equations, the long time-averaged one-point correlations 
of the large- and small-scale motions in the mean equation have to be provided. Here, the 
large-scale motions are calculated by solving their governing equations. These equations 
can be obtained by subtracting the mean equations from the short time-averaged Navier- 
Stokes equations: 

= o M 


dx i 


du{ 

dt 


+ 





1 d 7 ui 
Re dxjdxj 


-w > -«> 


dp 


dx, 


(6.6) 


Gaster et al. (1985) argued that the nonlinear terms can be neglected and showed that 
the local characteristics of the large-scale fluctuations in mixing layers can be described 
surprisingly well by linear theory. In fact, some weakly nonlinear solutions in hydrodynamic 
stability, for example, Stuart (1958), assume that the local shapes of the finite amplitude 
disturbances are those obtained by linear theory. Therefore, it is assumed here that locally 
the large-scale turbulent structures may be described by linear analysis and that their 
behavior is only weakly nonlinear. Hence, all the nonlinear interaction terms are neglected 
in the present formulation. 

The next assumption is that molecular viscous effects are unimportant. Davis and 
Moore (1985), in a numerical study of plane and axisymmetric mixing layers, found that 
the effect of decreasing the Reynolds number is to smear the vorticities without altering 
the dynamics of the large-scale structures. This phenomenon can also be observed in the 
experimental results of Brown and Roshko (1974) and Konrad (1976). In a mixing layer 
calculation, Tam and Chen (1979) showed that for a Reynolds number over 300, based 
on the local width, the unstable waves are not affected by the Reynolds number. In fact, 
increasing the Reynolds number produces more small-scale structures without significantly 
altering the dynamics of the structure of large scales. "Viscous effects are thus of minor 
importance in the development of the large-scale structure and are, therefore, neglected in 
the present approach. Computationally, this approximation also means a huge saving m 
computer time. 

The equations for the large-scale fluctuations can be simplified further by assuming 
that the mean flow is locally parallel. For free shear flows like mixing layers, wakes and jets, 
mean flows diverge slowly and this renders the locally parallel flow assumption applicable. 
For two-dimensional problems, the equations governing the large-scale unsteady turbulent 
motions, after introducing the above three physical assumptions, become 


du dv 
dx dy 


(7. a) 


4 


du 

~di 


dU 

~ V— 

= 

dp 

(7.6) 

v &y 

dx 


dx 


) 

: ~ u 

dv 

_dp 


(7.c) 

dx 

dy 



de not 

only the information that 

is necessary to close 


dv 

Tt 


X c^uauviiJ. , j - 1 U 1 t 

the system of equations for the mean flow, but also the behavior of the large-scale turbulent 

m0tl Si 1 n C e the coefficients in equation (7) are functions of y only, a simple separable form 
of solutions may be assumed, 

{u, V, p} = A{x) 9?{[u(y),0(y),p(y)] exp |*'(ai - w*)]) ( 8 ) 


where bold face quantities indicate a complex solution whose real parts lead to real quan- 
tities and SR denotes the real part of a complex number, a is a complex wavenumber an 
w represents the frequency. These are the normal mode representations generally used m 
hydrodynamic stability. The amplitude function A(x) appears as a parameter o 
calculations and is determined by the large-scale kinetic energy equation. This wea y 
nonlinear approach is usually referred to as a "wave envelop” method. In free shear flows, 
large-scale structures occur randomly in space and time and propagate ownstream m 
the form of quasi-periodic travelling waves. The wave-like solutions, equation ( ), us 
represent physical phenomena as well. The resulting equations obtained by substituting 
these expressions into the linearized inviscid equations for the large-scale disturbances can 
be reduced further to a second order ordinary differential equation in terms o v . 


{(“tf - ")(^i - = 0 (9) 

Equation (9) is the incompressible Rayleigh equation, Rayleigh (1880). 

For a given mean velocity profile and appropriate homogeneous boundary conditions, 
equation (9) forms an eigenvalue problem. In the present analysis, the complex wavenum- 
ber, a, of the disturbance is the eigenvalue. Note that the wavenumber appears non mear y 
in the problem. In addition to a traditional shooting procedure, three global approximation 
methods are applied in the solution of equation. These include the Chebyshev spectral 
approximation, a psuedospectral Chebyshev collocation method and a finite difference 

method. . . 

The amplitude function A[x) introduced in equation (8) along with the conventional 

normal mode representation determines the amplitudes of the coherent fluctuations. In th 
present analysis, the amplitude is determined by the kinetic energy equations for the large- 
scale fluctuations. Therefore, instead of growing exponentially, exp(-2a>), as would be 
predicted by linear theory, the development of the amplitude function is determined by t e 
conservation of the kinetic energy of the large structures. Equations for the kinetic energy 
of the large-scale fluctuations can be obtained by multiplying the momentum equations 


o 



for u„ equation (6.b), by u, and long time-averaging the resulting equation. The energy 
equation for the large-scale fluctuations can be written as: 


r . dk 

Oy— = ~u,Uj 
dxj 


,“£ - + 
J 5xy dzj 




— — — < ti 7 ti 7 >) + viscous terms (10) 

dzj ' ' } 

where k = iuTuy. Since the large-scale structures are inviscid in nature, the terms invo lving 
viscosity are negligible. Production of the large structure kinetic energy is positive if -u.uy 
and dUi/dxj are of the same sign, negative otherwise. Regions of “negative production 
associated with the large-scale coherent structures have been observed, Komori and Udea 
(1985) and Hussain and Zaman (1985). Conventional eddy viscosity models fail to predict 
this phenomenon. Since the basic assumption of these models is that momentum exchanges 
axe proportional to local mean flow gradients. That is, eddy viscosity models predict co- 
gradient momentum transport, or positive production of turbulent kinetic energy by their 
nature. In fact, the dominant structures in free shear flows are of large scale. The large- 
scale structures transport fluid elements across the flow unmixed and this is not directly 
related to local mean flow gradients. Energy is subsequently extracted from the large-scale 
motion and dissipated at the high frequency end of the wave number spectrum. The terms 
containing the residual stress tensor, - < u(u' >, describe the draining of energy from 
the waves. Very little information, experimentally or numerically, is currently availa e 
regarding these stresses. 

In the following section, three models are derived from the weakly nonlinear theory 
described above. The weakly nonlinear theory is formulated for turbulent free shear flows 
in general. As a test of the theory, these three models are derived for turbulent free mixing 
layers. The large-scale structures observed in mixing layers present most of the features of 
organized structures observed in turbulent free shear flows. A sketch of a turbulent mixing 
layer is shown in Figure 1. 


II. b Model I 

Free shear layers that possess inflectional-point instability, are inherently unstable. 
As the flow' develops, turbulence and/or background random noise provide perturbations 
necessary to excite these unstable waves and promote initial vortex formation. Therefore, 
the large-scale structures are made up of all modes residing in the flow. A complete simu- 
lation of the large-scale turbulence spectrum would require the inclusion of a broad range 
of frequency and spanw'ise wavenumber components. This was accomplished in the local 
solution of Tam and Chen (1979) and the integral model of Morris and Giridharan (1990). 
Figure 2 shows the unnormalized Reynolds shear stress distributions from Liou (1986) for 
waves of various frequencies. The velocity profile of the basic flow is a hyperbolic tangent 
function and the most unstable frequency for this flow is 0.2. The stress distributions are 
normalized by their peak amplitudes. It can be seen that for a wide range of frequencies 
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around the least stable mode, the distributions of flow quantities due to the large-scale 
fluctuation are similar. However, the mode that interacts most strongly with the mean 
flow is the one that has the largest spatial growth rate: the most unstable wave. Thus, or 
efficiency, instead of including all the unstable waves, it is assumed here that the o 
least stable modes are most effective in extracting energy from the mean flow. The most 
unstable waves are then used to describe the overall behavior of the coherent, large-sea e 
motions. In other words, the large-scale motions are described by the locally most amp i- 
fying disturbances in the flow. The method used here to locate the most unstable modes 

can be found in Liou (1986). _ 

The Rayleigh equation, equation (9), governs the local distribution of the large-sca e 

velocity fluctuation in the y direction. The equation is solved locally at each streamwise 
location. The amplitude function A{x) is determined from the energy equation for the 
large-scale structures, equation (10). 

The equation for the amplitude function is obtained by substituting the wave form 
expressions, equation (8), for the fluctuation components into the wave kinetic energj 
equation and integrating across the flow. The resulting equation is 

iG ' A ' = G,A> + C 3 (U) 


where 


Gi = - f + [C7(tt&* +t»v*) + 2.0 9?(up‘)J dy 
4 J —CO 

— 3?(ut> )J dy 


G 2 - -- 


where an asterisk denotes the complex conjugate. The terms containing the residual 
stresses in G 3 are responsible for draining energy from the waves. G 3 is of crucial im- 
portance in determining the wave amplitude and has to be considered carefully. However, 
very little experimental or theoretical results are available regarding these stresses. In the 
following analysis, we assume that the energy transferred out of the large-scale fluctuations 
is proportional to 

li! (12) 

l 

where u and l are the characteristic velocity and length scales of the large-scale motions. 
This model assumes that the turbulence is in an equilibrium state for the small-scale 
fluctuations. That is, the rate at which energy is transferred from the large scales is equal 
to the rate at which the energy is dissipated by viscosity at the small scales. The net effect 
of these terms may thus be modeled by 

Ci <»> 

where Cj is a model constant. The energy drain integral, G 3 , becomes. 

dy 


Ci 

Gs = --4 


r 

[,... , - «. . \ 3/2* 

(uu ■+■ vt> ) 

J — CO 



(14) 
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The amplitude of the wave structures can be calculated using an explicit fourth order 

Runge-Kutta method or the implicit Euler method. 

To solve the mean flow equations, the small-scale Reynolds shear stress —u'v' has to 
be modeled. Since it is mainly the motion at the large scale that is considered here, a 
simple zero-equation model is employed. That is, 


, ,dU u dU 

-U'V' = c, /- Igj-K aj- 


(15) 


The model introduces an additional parameter, C 2 . 

Thus, Model I contains two model constants: C\ used in the modeling of the energy 
transfer term in the amplitude equation for the wave structures and C 2 used in the simple 
eddy viscosity model for the small-scale Reynolds stress, u v . 


II.c Model II 

The large-scale structures in turbulent mixing layers are dynamically active and dom- 
inant. Thus the development of the flow is mostly affected by turbulent motions of large 
scale. Consequently, in this approach, only fluctuations at the large scale are included. 
Thus, there is no direct interaction between the small-scale structures and the develop- 
ment of the mean flow. This is also suggested by the analyses of Tam and Chen (1979) 
and Morris and Giridharan (1989). The characteristics of the locally most unstable mode 
is still considered to be representative of that of the large-scale structures. The energj 
gained by the large-scale structures from the mean flow is prevented from becoming un- 
bounded through the energy transfer terms in their energy equation. Thus, the influence 
of the small-scale turbulence on the mean flow is indirectly through the energy balance of 
the large-scale motion. In this model, the equations for the mean flow and the large-scale 
structures are the same as those in Model I. The model contains only one model constant 
, Ci, that describes the transfer of energy from the large to the small scales. Nevertheless, 
both Model I and II predict the average behavior of the shear layers. 

II. d Model III 

Model III simulates the time-dependent motion, at the large scale, associated with the 
passage of a train of large-scale structures. Experimental observations suggest that, even 
if initially there exists a continuous spectrum of infinitesimal disturbances upstream of the 
splitter plate, a disturbance emerges dominating over other neighboring perturbations in 
the early stages of the flow development. As the flow evolves; however, there is a contin- 
uous shift of the dominant component toward lower and lower frequencies. In fact, the 
growth of an initially small periodic disturbance is often followed by the development of 
subharmonics. In numerical simulations, however, the initial conditions can be conceived 
in a much simplier way. Instead of monitoring the disturbances in the initial continu- 
ous spectrum, a hierarchy of disturbances made up of the initially most unstable mode, 
according to linear theory, and its subharmonics may be chosen. This reflects the subhar- 
monic evolution modeP proposed by Ho and Huang (1981). Thus, the unsteady turbulent 
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large-scale fluctuations are described by the superposition of the instability waves in this 
hierarchy. This enables the time-dependent flow field at the large scale to be simulated. 
The solution methods for these equations are the same as in Model 1 and II. The equation 
for the amplitude function, however, has to be modified. Firstly, it is assumed that inter- 
actions between harmonics are negligible as there should be sufficient phase jitter in the 
unexcited shear layer. In addition, the details of the process of energy transfer from the 
large to the small scales is not modeled explicitly. At each axial location where a given 
instability wave saturates, or begins to decay, the energy is immediately remove rom 
the system. Consequently, there is no need to specify either a constant associated with 
the energy transfer process or the effects of the interaction between the small-scale motion 
and the mean flow. It should be noted that the energy equation is only solved for each 
instability wave during its period of growth. For amplifying waves the comparison between 
linear theory and experiments by Caster et al. (1985) showed that an inviscid analysis 
is adequate. Thus, the interaction between the large and small scales will be neglected 
during the growing or unstable region for each instability wave. 

The solution for the turbulent fluctuations is then fed back into the iteration process to 
get the corrected mean flow solutions. A visualization of the unsteady flow’ field predicted 
by Model III is made by means of streaklines. The streaklines are produced by injecting 
passive marker particles at the initial location, x = xo at various points across the shear 
layer. The positions of these particles at subsequent times can be calculated using the 
equations: 

= U [ *(t) , y{t) } + u [ x(t) , y(f) , t } ( 16 ) 

at 

and 

-j-y(i) = V | x(t) , y{t) ) ) + v [ x{t) , y{t) , t ], ( X/ ) 

at 

with 

x(0) = x 0 , yjfe(O) = yo , k = 

Particles thus move at each time step according to the local instantaneous velocity field. 

This concludes the description of the wave models. The dominant large-scale coher- 
ent turbulent structures in turbulent free shear flows are modeled in a weakly nonlinear 
manner. Three models are derived to simulate the development of turbulent free mixing 
layers. These models connect the development of the mean flow field with the oynamics 
of the large-scale turbulent fluctuations. The equations governing the mean flow field and 
the unsteady large-scale turbulent motions form a closed system of equations. 


III. Numerical Procedure 

The boundary-layer approximation renders the system of equations governing the mean 
flow parabolic. A fourth order Keller-Box scheme is applied to solve the resulting equa- 
tions. The equation for the instability wave, which is the Rayleigh equation in the present 
formulation, has been solved using various methods, including a traditional shooting, two 
spectral and a finite difference methods. For spatial instability, the system of equations 
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generated by the global approximations of the Rayleigh equation forms an eigenvalue 
problem which is nonlinear in its parameter, the wave number. It may be solved using the 
Linear Companion Matrix method or a method based on matrix factorization, Bridges and 
Morris (1984). Details of the various solution schemes can be found in Liou (1986,1990) 
The Rayleigh equation and the equations for the mean flow are solved iteratively at each 
streamwise location. The convergence criterion for the iterations is 



; 


where is a small number and M is the total number of grid points at each streamwise 
location. The amplitudes of the waves are calculated explicitly using a fourth order Runge— 
Kutta method applied to the wave energy equation (ll) 


IV. RESULTS AND DISCUSSIONS 


The models have been tested in an incompressible free mixing layer. A hyperbolic 
tangent distribution is taken as the shape of the initial streamwise mean velocity, U(zq, y), 
i.e., 

U ( x 0 , V ) = \ ( 1 + tanh ( 30 y )). (19) 

The initial cross-stream mean velocity, V(xo,y)j can be se ^ a sma ^ number or zero. 
The boundary conditions for the mean flow are 

U ( x , y u (x) ) = l.o, U ( x , y/(x) ) = 0.0, V ( x , y u (x) ) = 0.0. (20) 

where y u (x) a nd yi{x) are the upper and the lower boundaries of the physical domain shown 
in Figure 3. As a test of the ability of the instability wave model to describe large-scale 
structures and the associated Reynolds stresses, the model was first applied in the self- 
similar region of the flow with a mean velocity profile from Patel (1973). Figure 4 shows the 
calculated and experimental Reynolds shear stress distributions. Calculated results using 
a traditional shooting method compare favorably with that using global approximations 
of various order. Note that all the calculated results have been normalized by the peak 
experimental value. The discrepancy at the low— speed side of the layer suggested that the 
momentum exchanges due to the small-scale turbulent motions might not be negligible 
in this region. It should be noted that this negative value of Reynolds shear stress disap- 
pears for small values of velocity in the lower stream. The structures obviously contribute 
negative shear stress at the low-speed edge of the flow. This counter— gradient transport 
of momentum gives negative energy production in this region. A similar phenomenon was 
observed experimentally by Komori and Ueda (1985) in the self-similar region of a jet. 
In fact, regions of negative shear stress can be easity observed if the large-scale struc- 
tures are excited artificially, for example, see Wygnanski, Oster and Fiedler (1979) This 
counter— gradient momentum transfer decelerates and subsequently reverses the flow on 
the low— speed side of the mixing layer as the shear layer develops. 


10 



IV. a Model I 


As noted above, Model I proposes that a contribution from the small-scale Reyno s 
stresses is required to describe the total turbulent forces that determine the developmen 
of the mean flow. The model introduces a new parameter, C 2 . Latigo (1979) argue a 
the turbulent shear stress contributed by the small-scale, incoherent motions is about a 
half of the total shear stress. An estimate of C 2 based on the value that is use in e 
classical eddy-viscosity models is then obtained. In addition, the force terms associated 
with the large-scale normal stresses in the mean momentum equations are also retail e , 
since they are found to be of the same order as the other Reynolds stress gradient terms on 
the low-speed edge of the shear layer. The normal stresses associated with the large scale 
structures can be calculated directly by the wave models and involve no further empiricism. 


In the numerical calculations, the local solution of the Rayleigh equation is found 
time-consuming. To accelerate the axial marching an adaptive grid has been evise 
grid size in the cross-stream direction in the transformed domain are fixed. The axial step 
sizes are chosen such that the convergence indices of the first iteration at a downstream 
station are greater than a fixed number e 2 . 



; 


(21) 


The grids are found to cluster in a region where there are large changes of flow properties 
for example, when the flow is developing initially. 


The initial wave amplitude represents the initial strength of the instability waves or large 
scale motions for which there are no quantitative experimental measurements An estimate 
of the initial wave amplitude can be made based on the initial energy flux of the tur u enc 
From a sequence of numerical experiments, however, it is found that flows with relatively 
strong initial amplitude saturate early. Subsequently the flow develops in a similar manner 
only the virtual origin of mixing is changed. The initial amplitudes for the cases presente 
in this paper are fixed at 5 x 1CT 2 . The corresponding initial turbulence intensity is about 
1%. The model constant C x of the energy transfer term in the wave kinetic energy equation 
is taken from a conventional Prandtl energy model, Launder, et. al. (1979). t is oun 
that its value has no significant influence on the results of the mixing layer calcu ations 
values of C\ and C 2 used here are 

Cl = 2.8 , c 2 = 0.08 ( 22 ) 


Again, in Model I, turbulent forces associated with the wave shear and normal stresses 
as well as the small-scale motions are considered. Figure 5 shows the axial forces acting 
on fluid elements across the layer at various axial stations. Negative or retar mg orces 
associated with the wave Reynolds shear stress appear near the zero speed side early m the 
developing stages of the layer and never change sign as the flow develops. . evert e ess, 
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the other two mechanisms, which are the gradients of the normal stresses of the turbulent 
large-scale structures and of the shear stress of the small-scale turbulent fluctuations, 
contribute positive driving forces. The net effect is the co-gradient momentum_transport 
near the zero velocity side of the layer. It can also be seen that the — (u — t ) x terms 
plays a significant role in the dynamics of the mixing layer development, especially near 

the zero velocity side of the layer. 

The growth of the layer measured by 6 is shown in Figure 6 6 is the distance between 
the points where the local mean velocity is 0.9 and 0.1 of the main stream mean velocity , 
U. The calculated rate of growth, dS/dx, agrees reasonably well with the value that is an 
average taken over various experiments in the self-similar region of the flow . 

Figure 7 shows the predicted axial mean velocities at a sequence of downstream sta- 
tions as functions of the self-similar coordinate, r?. The agreement between the predictions 
and experimental data is good except near the low speed side of the shear layer. The 
laver has not reached self-similar state at i = 0.72 and local velocity profile differs from 
Patel's self-similar velocity profile. Note that for a free mixing layer, the accuracy of the 
measured mean flow data in the low speed region is poor due to the rapid variations m 
the instantaneous flow direction. Any agreement between results on the zero speed edge 

is likely to be fortuitous. 

As can be seen from Figure 8, which shows the shear stress distributions across the 
mixing layer at various stations, the sum of the shear stresses from the large sea e an 
small-scale motions agrees well with experimental data. The experimental measurements 
are the long time-averaged correlations of the turbulent fluctuations. Model t us not 
only predicts the mean velocity profiles but also appears to model correctly the turbulent 
momentum transport in the layer. The latter is usually achieved only by using lg er 
order moment closure models, which includes a large number of model constants. 

The amplitude of the large-scale fluctuations is plotted in Figure 9. The large-sea e 
structures extract energy from the mean flow and grow as the flow develops. However, 
energy is also being transferred to the smaller scales and subsequently dissipated by viscos- 
ity. The final equilibrium of the large-scale motions amplitude is reached when the energy 
gained from the mean flow balances the energy lost to the small scales. 

FV.b Model II 

Model I assumes that the large-scale and the small-scale turbulent fluctuations play 
a direct role in the momentum transport process in the mixing region. Therefore, m 
Model II only the fluctuations at the large scale are included. This eliminates the need to 
specify a model, equation (15), to describe the momentum exchanges due to the small-scale 
fluctuations. Model II thus involves only one model constant for the energy dissipation 
model in the kinetic energy equation for the large-scale turbulent fluctuations. It can 
be seen from Figure 10 that the forces associated with the large-scale normal stresses 
are apparently able to counter-balance the decelerating effects of the wave shear stress 
gradients. Thus, a prediction of the development of the layer is possible considering only 
the dynamics associated with the large-scale turbulent fluctuations in the layer. 
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The predicted mean velocity profiles are presented in Figure 11. It shovs that the 
mean flow can be satisfactorily predicted by modeling only the dominant larg^e.tm- 
tures Figure 12 compares the calculated axial mean velocity profiles using 
II Note that the predicted development of the mixing layer has reached equilibrium state 
both x = 6,9 L, Mode, 1 and x = 7.26 os.ng Mode, II. Thus, Figure ,2 
the mean velocity profiles at the equilibrium state of a mixing layer usi g 
turbulence closure schemes. The agreement between these predictions is more than satis- 
factory. Both of the predicted mean velocity profiles deviate from Patel s data near the 
low speed side of the laver. However, as noted previously, the accuracy of the measur ^ 
data may be suspected in that region due to the rapid variations in the ^tantaneous w 
direction. This phenomenon, in fact, is predicted in the application of Mo e 
turbulent mixing layer. 

As is shown in Figure 13, the predicted shear stress distributions do not milch e 
total shear stress distributions measured by Patel (1973). However, as noted above this 
difference does not necessarily mean that the small-scale stresses should be me u e • 
must be remembered that the present model simulates the entire large-scale spectrum with 
a single frequency wave that is locally most unstable. Tam and Chen (1979), m their local 
model, included a broad range of instability waves and found good agreement with exper- 
iments without the inclusion of contributions from the small scales. In fact, it is sftown 
here that the time average characteristics of a turbulent mixing layer can be predicted sa ■- 
isfactorily using only the most unstable waves, provided that all the momentum transpor 

mechanisms associated with the wave are included. 

The evolution of the large-scale amplitude using Model II follows a similar behavior 
to that using Model I and is shown in Figure 14. Once again an equilibrium condition is 
reached where the rate of energy transfer from the mean flow to the large-sca e struc ures 
balances the rate at which energy is lost by the structures to small sea es or even ua 
dissipation. The little kink near x = 6.0 is due to the fact that at this region there are 
relatively large changes of marching step sizes required by the adaptive grid genera ion 
scheme. Since the marching step sizes are selected based on the globa variation o 
mean flow, this kink has little effect on the predictions of the flow development. 

Figure 15 gives the growth of the mixing layer in terms of momentum thickness, , 
predicted by Model II. In the present analysis, the momentum thickness is defined as 


*= r 9 ( l-,)dy ( 23 ) 

Jlyi 

6 0 is the initial momentum thickness and the straight line represents the rate of growth of 
the layer, which is an average over experimental data at the equilibrium state o mixin 0 
layers. In the earlv stages of the development of the mixing layer, the large sea e str 
in the flow are relatively weak and momentum exchanges are mainly due to the eftect o 
molecular viscosity. The strengthening of the large-scale structures increases drastically' 
the mass and momentum exchanges across the layer and, consequently-, the wi o 
mixing region. This reflects the same phenomenon predicted using Model I. 
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Figure 16 shows the Strouhal number, Ste , of the large-scale structures based on the 
local momentum thickness. That is, 



u>6 

2nU 


(24) 


where / and w denote the frequency and the angular frequency. After the mixing layer 
reaches an equilibrium state, the Strouhal numbers of the large-scale distur- bances ap- 
proach a constant (« 0.012). The Strouhal number of the large-scale structures base on 
the average passage frequency and the local momentum thickness, in an unforce , mi la ) 
fully turbulent mixing layer is «0.024, Hussain and Zaman (1985). In a spatially devel- 
oping mixing laver, unstable waves or large-scale structures are continuously amplified as 
they propagate downstream. The amplification of the unstable waves continues until they 
become neutral. Thus, a wave at its neutrally stable stage reaches maximum amplitu e 
and dominates over waves of other frequencies. Consequently, the detected average pas- 
sage freouency of the large-scale structures is associated with that of the locally neutral 
mode, which is about two times of that of the locally most unstable mode. The present 
calculation, which predicts that the Strouhal numbers of the locally most unstable mooe 

reflects this phenomenon. 


IV. c Model III 

Model III simulates the time-dependent turbulent motion associated with the passage 
of a train of large-scale turbulent structures. The large-scale turbulent structures are 
represented by a superposition of hydrodynamic instability waves. As the flow oevelops 
axially, these hydrodynamic waves become damped because of the growth of the shear 
layer. Since it is assumed that energy associated with a given wave is removed immediately 
after it becomes neutral, there is no need to obtain damped inviscid solutions by analytic 
continuation in the complex plane, Tam and Morris (1980). 

The initial wave amplitudes of this calculation are 

Aoy = 10" 2 , J = 1.--6 ( 25 ) 

The initial mean velocity profiles and the boundary conditions are the same as those used m 
the previous calculations. In the preliminary calculations, it was found that an abnorma ltv 
in the mean velocity distributions appeared near the critical points of saturating v ‘^' s - 
Also, most of the shear layer growth occurred on the low speed side of the layer. is 
gives a non-monotonic velocity distribution near the critical layers of saturating waves and 
another inflection point appears. Saturating waves thus have to be removed before they 
become neutral during the axial marching. Wygnanski and Petersen (198/) suggested tha 
this abnormality is due to nonlinearity. Composite expansion techniques have been applied 
to investigate the effect of critical-layer nonlinearity, for example, see Goldstein and Leib 
(1988) and Goldstein and Hultgren (1988). Another approach to resolve this issue is to 
include viscous effects: that is to solve the Orr-Sommerfeld equation. Since the present 
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investigation is directed toward developing simple turbulence models, instead of including 
other computationally expensive approaches, the effects of the critical point is accounte 
for by incorporating a small amount of eddy viscosity in the analysis of the mean ov . 
The additional mixing at the fine scale is diffusive and able to smooth out humps m e 
flow. In the present calculation, the extra mixedness provided accounts for about 10/c ot 
the amount of turbulent momentum exchange that is suggested by conventional models. 

With this modification the mean velocity distributions predicted at several dov m 
stream stations are shown in Figure 17. There are six waves in the hierarchy in this 
calculation. Since waves are removed successively during the axial marching, t e num er 
of waves included depends on the distance the calculation is to be carried downstream. 
There are some small differences between the calculated results and Patel's measurements. 

It should be noted that Model III simulates the development of the mixing layer asso- 
ciated with the realization of a single event, which is the passage of a train of well-defined 
large-scale structures. In physical experiments, such as that reported in Patel (19 73 )> ex- 
ternal disturbances or noise from various sources may modify the initial conditions o 
flow, which cause events that affect the subsequent development of the layer. There ore, 
the differences between predictions using Model III and experimental data will vanish if 
these randomizing events are taken into account. On the other hand, it mig t e ex 
pected that Model III should resemble the behavior of that of an externally excited mixing 
layer, in which the layer is excited at particular frequency and cleaner flow pictures can 

be observed. 

Figure 18 shows the development of the wave amplitudes. The additional small-sca e 
mixing increases the initial growth of the layer so that the fundamental mode is remove 
at a lower amplitude than its subharmonics before its amplitude reaches equilibrium level. 

The axial width of the layer is shown in Figure 19 and is compared with the prediction 
using Model I. As was noted earlier, the presence of this stepwise evolution is characteristic 
of excited flow's and would be smoothed out if many waves with slightly different amplitudes 
and frequencies were included. However, the global evolution of the width of the layer 
agrees closely with that using Model I. Figure 20 compares the predicted evolution of 
the momentum thickness of the shear layer, using Model II and Model III. Note that 
the initial momentum thickness for the cases using Model I, II and III are the same 
The only differences are the values of the initial amplitudes. The case using Model III 
assumes stronger initial large-scale structures than that used in Model II. As the flo 
develops downstream, Model III predicts a greater amount of large-scale mixing of mass 
and momentum than Model II. Consequently, the predicted grow'th of the mixing region 
using Model III is faster than that using Model II. However, the predicted rates of growth, 
d0/dx, using these two models are virtually the same. Since both of these models provide 
predictions by modeling the intrinsic characteristic structures in the free mixing layer, it 
is not surprising to find some family resemblance between the results predicted by these 
three models. 

Figure 21 shows the unsteady velocity profiles in the axial direction at x/<5 0 =30 before, 
during and after the passage of a large-scale structure. At this location, the dominating 
mode has the period of approximately 4. Reverse flow occurs at the low speed side o 


the layer during the passing of a large-scale structure. At the same time, the streamwise 
velocity increases instantaneously at the high speed side of the layer. These instantaneous 
velocity variations show the passage of a clock-wise rotating structure. Note that the 
instantaneous reverse flow at the low speed side of the layer gives rise to instantaneous 
changes of flow angles and makes velocity measurements extremely difficult in that region. 
The visualization of the flow can also be assisted by streakline plots such as those shown in 
Figure 22. The roll-up of vortices into larger vortex-like structures can be observed clearly. 
The initial roll-up is dominated by the fundamental mode. As time progresses* the initial 
structures convect downstream and roll around each other. These regions of concentrated 
vorticity then form a single large structure. As the passive particles travel downstream* 
their motion becomes dominated by lower subharmonics. Vortex-like structures of increas- 
ing scale are formed. Subsequently, the rolling process between two adjacent structures 
repeats as the flow develops further downstream. Careful examination of the figures also 
shows how the structures are convecting downstream as they form and roll. Large tongues 
of unmixed fluid are swept across the layer and reach the opposite side of the layer as ob- 
served by Brown and Roshko (1974). The engulfed fluid elements from the two sides of the 
layer mix and are drawn into the leading and trailing vortices when passing through the 
high-strain braid region between the vortices. This provides the environment for further 
fine-scale mixing. 

Figure 23 and Figure 24 show the flow pictures frozen at t = 6.5 and 5.0, respectively. 
In these cases, harmonic waves in the wave hierarchy are in phase. The distribution of 
momentum thickness in the streamwise direction at t = 6.5 is shown in Figure 23. b. At 
t = 6.5, three full-grown large-scale structures centered at roughly, x/6q = 10,18 and 40 
can be observed. These large-scale structures are essentially vortices rotating in a clock- 
wise manner. The mean and transient velocity profiles at x/6o=lS and 40 are shown 
in Figure 23. c. The turbulent large-scale structures, which appear as clock-wise rotating 
vortices, contribute velocity excess/deficit on the high/low speed side of the layer* relative 
to their respective mean velocity distributions. Therefore, the instantaneous momentum 
thickness of the flow in the region occupied by the fully-grown large-scale structures is 
smaller than the mean value. This can be clearly seen in Figure 22. b. Note that for a 
vortex sheet with a velocity profile given by 
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i, v > o 

o y < 0, 


(26) 


the momentum thickness is zero. The mean (or instantaneous) momentum thickness at 
any axial station is obtained by substituting u (or < u > ) into the g in equation (26). In 
Figure 23. b* three dips can be observed clearly in the instantaneous momentum thickness 
distribution. The positions of the dips correspond to the centers of the three fully-grown 
structures. On the other hand, the instantaneous velocity profile at x/6q = 27 shows 
velocity deficit/excess on the high/low speed side of the layer. At i — 6.5, the high- 
strain braid region between two structures passes through z/6o = 27. Therefore it is not 
surprising to observe that the instantaneous momentum thickness is greater than the mean 
momentum thickness at x/6q — 27. If there is a very strong velocity excess/deficit on the 
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high/low speed side, a shear flow may have negative momentum thickness. The case s 
in Figure 24 reflects the situation. At t = 5, two vortices are merging at x/6 0 - ^ ine 
two merging vortices deform as they roll around each other. The instantaneous ve oci J 
profile at the center of the merging is shown in Figure 24.c. It shows clearly t e oot P nn ® 
of these two deformed structures. This results in high velocity excess/deficit on e lg 
low speed side of the layer. The local instantaneous momentum thickness becomes nega ive 
and reaches a minimum at xj6 0 = 30 before picking back up downstream. The following 
increase in momentum thickness in the axial direction is the effect of further mixing o 
low and high speed fluid. Similar to the case at t = 6.5 and x/6o—27, the instantane 
momentum thickness at x/6 0 = 30 and 50 are much higher than the respective mean values. 
As mentioned earlier in the analysis, harmonic waves have to be cut off before t e> satura e 
to avoid the problem associated with nonlinearity. Therefore, the rather abrupt 'aria ion 
in the instantaneous momentum thickness distributions in Figure 23. b and '| ure • 
can be observed at the axial locations where the waves are cut off. The second an t e 
third harmonics are cut off at x/6q = 18 and 30, respectively. The dotted lines in lg 
23 b and 24.b show the projected distributions of the instantaneous momentum thickness 
should the harmonic waves be carried through their neutral points in the calculations. 

V. SUMMARY 

Three models based on a quasi-linear theory, that describes the dynamics of the domi- 
nant large-scale structures in a free mixing layer have been presented. The closure schemes 
incorporating the models are able predict the development of the turbulent free mixing 
layer accurately, even though they contain some assumptions and simplifications I e 
predicted averaged properties of the incompressible turbulent mixing layer agree well wi 
experiments. The transient turbulent motions at the large scale in the layer mappe ou 
using Model III possess many features that are apparent in flow visualization experiments, 
such as the convective nature of the large-scale structures, the large-scale transport of un- 
mixed fluid elements and the roll-up of vortices. The models involve less empiricism t an 
most conventional models. Since large-scale coherent structures appear also in shear flows 
of other geometries, the closure schemes presented here should be applicable to those cases 
as well. It is hoped that these models, which originate from observed physical phenomena, 
will provide efficient tools to model other free shear flows. 
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Chapter 4 


Conformal Grid Generation Using Wegmann’s Methods 


Roy S. Baty and P. J. Morris 






ABSTRACT 


Conformal coordinate transformations are used to map simple computational 
domains onto arbitrary simply and doubly connected regions with smooth bound- 
aries. Efficient schemes involving the solution of the inverse boundary correspon- 
dence function problems associated with mapping the unit disc or circular annulus 
onto simply or doubly connected domains respectively are employed. The numer- 
ical implementation of these schemes is emphasized. Examples are generated for 
regions with elliptic inner and outer boundaries. Additional examples are used to 
demonstrate the accuracy and convergence of the schemes and their practical lim- 
itations. The techniques are found to converge well if holomorphic functions are 
used to describe the boundaries. The use of preconditioning maps is also discussed. 
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1.0 INTRODUCTION 


Many problems of practical interest in engineering involve the solution of dif- 
ferential problems in complex geometries. In such problems the boundaries may 
not conform to coordinate lines in an orthogonal coordinate system. Alternatively, 
the coefficients in the differential problem may not be constant along coordinate 
lines making a separable form of solution impossible. Several numerical techniques 
are available to overcome these difficulties. Finite element or boundary element 
methods may be used. If a finite difference or spectral approximation is sought the 
physical domain must be transformed into a simple computational domain. Once 
again, several alternative approaches exist. The use of conformal maps is very 
desirable. Such maps simplify the governing differential equations in the mapped 
regions since the metric tensors are diagonal. However, these maps are difficult, 
and sometimes inefficient, to generate computationally and often lead to ill-posed 
numerical problems. 

Recently, Wegmann [1,2] developed a very efficient scheme to determine the 
conformal map from a standard computational domain onto an arbitrary simply 
connected region in the plane. This scheme solves the inverse boundary correspon- 
dence problem associated with mapping the unit disc onto a region with a smooth 
boundary. Wegmann [3] also extended this technique to determine the boundary 
values for the transformation mapping the circular annulus onto a doubly connected 
region with smooth boundaries. In both cases, the entire conformal map may be 
generated from the solution of the boundary correspondance problem using the 
Cauchy Integral Theorem. 

The present study of conformal maps is motivated by the authors’ interest in 
the spatial stability of jets of arbitrary cross section. The characteristics of the 
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instability waves play an important role in the jet mixing process and the radiation 
of noise by high speed jets. The growth of these instabilities is governed by the 
linearized, inviscid equations of motion. These may be simplified to a single linear, 
elliptic, second- order partial differential equation for the instability wave pressure. 

In general, the coefficients in this equation vary arbitrarily in the cross section 
plane and no separable solutions may be found. Thus a numerical solution must be 
obtained in the regions of non- constant coefficients. 

There are two characteristic regions in a jet flow. The first is the annular mixing 
region surrounding the potential core of the jet and the second is the developed jet 
flow region downstream of the potential core. Once a map is generated for either of 
these doubly or simply connected regions onto a simple computational domain, the 
homogeneous boundary value problem for the instability wave pressure fluctuations 
may be solved. Such a solution, using a hybrid pseudospectral and finite difference 
algorithm is described by Baty and Morris (4). 

In this paper we apply Wegmann’s techniques to compute the conformal co- 
ordinate transformations for simply and doubly connected regions. Examples are 
given for regions appropriate for the study of jets issuing from elliptic nozzles. Sim- 
ply and doubly connected elliptic regions of aspect ratios 2 and 3 are considered. 
The conformal maps for such regions are difficult to compute and provide a good 
numerical test case. The numerical implementation of the mapping techniques is 
emphasized. The elliptic and additional maps are generated to establish the con- 
vergence and accuracy of the techniques. In addition, the practical limitations of 
these techniques and some methods to overcome these difficulties are presented. 

In the next section some preliminary mathematics is developed. Sections 3 
and 4 describe the simply and doubly connected Wegmann methods. Details of the 
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numerical implementation of the techniques and several examples are given. Some 
practical limitations of the techniques and a discussion of convergence and efficiency 
is given in Section 5. 

2.0 PRELIMINARY MATHEMATICS 

The regions of interest in the present study are the annular and disc-like regions. 
Therefore, the canonical regions are the unit annulus (with some inner radius 0 < 
li < 1 ) and the unit disc. Let C denote either canonical region, let P denote a 
given physical region, and let W denote the conformal mapping satisfying: 

W:C-*P 2.1 

Consider a rectangular region, P, defined by 

£:= [(y\y 2 ) eR’Ia^y 1 <b,0<y 2 < c] 2.2 

where a, 6, and c are finite real numbers and where c > 0. The exponential mapping, 
exp(iz), then carries R onto an annular region, C : 

exp : R — > C 2.3 

Composing the mappings given by 2.1 and 2.3 yields a map from a rectangular 
domain onto the region of interest: 

W o exp : R — *■ P 2.4 

The exponential map carries an infinite strip onto the unit disc. Thus, any finite 
rectangle will be carried onto an annulus. 

The following discussion considers the analysis required for simply connected 
geometries. The analysis needed for the doubly connected geometries may be gen- 
eralized from this case. 
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Let V be the unit disc. Let i be a simply connected region with a smooth 
closed boundary parameterized by, 

z(0 on 0 < t < 0. - 2 - 5 

Here, z is assumed to be a smooth regular curve. By the Riemann mapping theorem 
there exists a unique function g such that, 

g : D -> i. 26 

It should be noted that the uniqueness of g follows from the imposition of normal- 
izing conditions, for example see [5]. For simply connected regions the following 
conditions are imposed: 

0(0) = o, ff'(O) > o. 2.7 

The goal of the simply connected Wegmann method is to determine the map g on 
the boundary of the unit disc. This requires determining the image of a point on 
the unit circle. Since <7(P) = <f , it follows that a point on the unit circle is carried 
to a point on the boundary of <£. Now, since the boundaries of D and £ are known 
smooth functions, only the angle of a point on the boundary of D or £ needs to be 
known in order to determine its location on the boundary. Therefore, if the angle 
0 of a point is given, the problem becomes to determine the angle r(6 ) of its image 
satisfying: 

»(«'•) = *(>■(*)) 28 


Any real function r{6 ) such that, 


is 2 7 T periodic and 2.8 is satisfied is called an inverse boundary correspondence 
function, IBCF. 

The Wegmann method solves for the function, r, so that equations 2.8-9 are 
satisfied. The technique assumes that a good guess, f, for r is known and then 
computes a small correction factor, t], such that: 

T = f + T) 2.10 

To simplify the construction of r, an approximation of t] is calculated. Substitution 
of 2.10 into 2.8 and linearization of the result yields: 

*(««) = z(m)+*'{no)M6) 2 - n 

where the prime denotes differentiation with respect to t. Equation 2.11 can then 
be recast as a function-theoretic boundary value problem allowing the explicit de- 
termination of rj. In the numerical solution of the IBCF, equation 2.11 is solved 
iteratively becoming: 

i(e‘ e ) = z{r k [6)) +z'M0))»? t (0) 2 - 12 

Therefore, at each step of the iteration the function r as well as the boundary of 
the smooth domain € are approximated. Once the IBCF, t, has been computed, 
2.8 may be used to construct the function, g , on the interior of the unit disc from 
the Cauchy Integral Theorem. 

3.0 THE SIMPLY CONNECTED WEGMANN METHOD 

Wegmann’s technique solves the inverse boundary correspondence problem, 
introduced in the last section, by iterative computation of the correction factor, rj , 
defined by 2.10. This method generates simultaneously two sequences of functions. 
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The first sequence of functions, g k , is analytic on the interior of the unit disc and 
satisfies the normalization conditions 2.7. The second sequence of functions, z t , 
maps continuously the unit circle onto the boundary of the physical region. If these 
sequences converge, they may be used to compute the boundary correspondence 

function, t. 

The basic strategy of the Wegmann method is to construct the real part of the 
desired conformal map, g kt on the boundary of the physical regions and then use a 
conjugate integral operator, K , to determine the imaginary part. Here K is defined 
by the equation, 

lfu(r) = ^P.V./ u(t)cot(^)dl 3.1 

where P.V. stands for the principal value of the integral. 

Equation 3.1 defines a general operator which generates the conjugate periodic 
function, v, of a real 2 tt periodic function, u. Furthermore, the operator K is easy 
to evaluate numerically. If the functions on which this operator acts are expressed 
as a Fourier series, K simply multiplies the coefficients of the Fourier series by ±i 
or by 0. The mathematical and numerical details are given in Henrici [6]. 

Wegmann’s method |3], constructs the real part of the conformal map by re- 
casting 2.12 as a Riemann-Hilbert problem given by: 



This Riemann-Hilbert problem has a unique solution, Henrici [5]. Furthermore, 
the solution of the Riemann-Hilbert problem may be used iteratively to compute 
the desired conformal map on the boundary of the physical region. The k + 1-th 
approximation of the map is given by : 

u , k (6) = K(0(T k (6))-8) 3.3 
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q k ( 6 ) = S[2(r*(0))exp(w k (0) - »>(r*( 0 )))] 

g k + l (e if ) = {iq k {6) ~Xk ~ Kq k {6))exp{i<t>{T k {0)) ~u k {6)) 

where <f> is the tangent angle of the curve defined by: 


«'(£) = |*'(0l ex P( l '^(0) 


Furthermore, the constant, Xk , is defined by: 

Xk = 9k cot a k 


where, 

* = hi Q ‘ [(]d( ' 3 ' 8 

and 

&k = J 0 3-9 

The constant X k insures that the normalizing conditions 2.7 are satisfied. 

The update of the correction factor, r) k , is then determined by 2.12 and 3.5, 
which yield, 


Vk{0) 


f z(r k (fl)) \ Xk + KOk (A) 

U'W 5 ))/ k'K( 0 ))i ex p( w *( 0 )) 


3.10 


Then using 3.8, the next iteration of the inverse boundary correspondence function 
becomes: 

T k + i{0) = Tk {0) + r)k{8)- 3,11 

Equations 3.3-11 constitute the simply connected Wegmann method. 
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3 . 1 : Numerical implementation of the Wegmann method 

The Wegmann method has been implemented numerically using FFT s to ap- 
proximate the Fourier series and the necessary periodic functions and their conju- 
gate integral operators. The FFT’s have been computed using the standard IMSL 
subroutines DFFTRF and DFFTRB. 

The present calculations have been performed on a Vax 11/780, in double 
precision, using the following steps: 

[1] An initial guess for the inverse boundary correspondence function, r„, is made. 

[ 2 ] The boundary functions 2 , z' and the argument of z ' , <f>, as defined by 3.6 are 
computed at the values of r k . 

[3] The integral equation 3.3 is then computed. This is accomplished using an FFT 
to approximate the integrand of 3.3. Then the conjugate integral operator is 
applied. Finally an inverse FFT is applied to obtain a discrete representation 
of the function u k . 

[ 4 ] Using the results of steps 2 and 3, g k is then computed from 3.4. 

[5] The conjugate integral operator is then applied to q k .This is accomplished using 
an FFT and its inverse is also obtained as outlined in step 3. 

[ 6 ] The normalization constant, defined by 3.7, is then computed, using 3.8 and 
3.9. 

[ 7 ] The final step is to use the results of steps 2 through 6 to compute r ] k , given 
by 3.10, and then update the inverse boundary correspondence function, 3.11. 

If a good initial guess is provided, the Wegmann method converges quadrati- 
cally. For the present application of grid generation, the scheme requires typically 
5 iterations. Therefore, the coupling of the FFT’s with the quadratic convergence 
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rate makes this scheme very efficient. 


3.2: Validation of the simply connected Wegmann method 

The code has been validated by computing the inverse boundary correspon- 
dence function for an inverted ellipse defined by: 

2 ( 5 ) = (1 - (1 - p 2 )cos 2 s)* exp(ts) 3.12 


for 0 < p < 1. In this case the boundary correspondence function may be deter- 
mined exactly. From Gaier (6], the solution of the inverse boundary correspondence 
function is given by: 

tan(0) = ptan(r) 3.13 


Table I shows the maximum value of the absolute errors for the inverse boundary 
correspondence function defined by setting p = 0.6 in 3.12. Here, N and k represent 
the number of points used to discretize the circle and the number of iterations 
respectively. These error results have the characteristic properties of the results 
found by Wegmann [l], for the inverted ellipse. Table I also illustrates the quadratic 
convergence rate of the scheme. 


The Wegmann scheme has also been used to compute the inverse boundary 
correspondence functions for regions bounded by ellipses. The boundary curves 
defining the ellipses are written in terms of holomorphic functions. These are given 
in Section 5. The sensitivity of the scheme to the functional form of description of 
the boundaries is also discussed in Section 5. Once the IBCF’s are determined, the 
conformal maps are computed using the Cauchy Integral Theorem, see Henrici [5], 
given by: 



3.14 
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Equation 3.14 is approximated by a series, using an FFT to compute the coefficients 
of the series. The series is then evaluated at the grid points in the computational 
domain by applying the Horner summation technique: Burden and Faires |7]. 

Figures 1 and 2 show the grids generated using the Wegmann method for simply 
connected elliptic regions of aspect ratio 2 and 3. Figure 3 shows the computational 
domain for these examples. Recall that this rectangle is mapped onto the elliptical 
regions by composing the exponential map with the Wegmann map. 

4.0 DOUBLY CONNECTED WEGMANN METHOD 

The doubly connected Wegmann method [3] is a generalization of the simply 
connected scheme. In this case, the canonical region for the conformal map is 
the unit annulus, with interior radius n . This scheme requires the solution of two 
inverse boundary correspondence problems computed simultaneously, one for each 
boundary of the physical region. Furthermore, the value of the interior radius of 
the canonical annulus must be determined. 

Since the Wegmann method for annular regions is a generalization of the simply 
connected method, the iterative steps are not outlined in the present study, but 
may be found in reference (3). It should be noted, however, that Wegmann tested 
two versions of the scheme. These versions differ in the iterative method used to 
compute the inner radius of the annulus, /i. All the numerical experiments presented 
in this section are based on Wegmann 's first scheme. The first method is applied 
since its convergence properties have been justified rigorously, while there is less 
mathematical justification for the convergence properties of the second scheme. 

The doubly connected scheme has been implemented numerically using FFT s 
to approximate the periodic functions and their corresponding conjugate integral 
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operators. As with the simply connected method, the calculation is initiated by 
providing a guess for the IBCF’s, r ll0 , and r 2 , 0 . The subscripts 1 and 2 denote 
the outer and inner inverse boundary correspondence functions respectively. With 
good initial guesses for the IBCF’s and the interior radius, the doubly connected 
Wegmann method converges quadratically. 

To validate the doubly connected Wegmann method, the inverse boundary 
correspondence functions have been computed for the annular region defined by the 
curves: 

21 ( 5 ) = exp(ts), 41 


and, 


z 7 {s) = A + B exp (is) 


4.2 


for 0 < A, B < 1. The conformal map carrying the unit annulus onto this region is 
given by 


m = 


jz + t 


4.3 


mz + n. 

where the constants j, £, m, and n may be determined from the coordinates of the 
center of the circles. Thus, the inverse boundary correspondence functions for this 
example may be computed easily from 4.3. 

Table 2 shows the maximum value of the absolute errors for the IBCF s defined 
for the non- -concentric region for A — 0.2 and B = 0.5. These results are again sim- 
ilar to the results of Wegmann |3] for this region. As in Table 1, N and k represent 
the number of points used to discretize the boundaries and the number of iterations 
respectively. The error results are presented in pairs: the first number corresponds 
to the error on the outer boundary, while the second number corresponds to the 
error on the inner boundary. 
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The Wegmann method has also been applied to compute the boundary corre- 
spondence functions for annular elliptic regions. The functions used to define the 
boundaries of the elliptic regions are given in Section 5. Once the IBCF’s are com- 
puted, the conformal maps are computed using the Cauchy Integral Theorem. For 
the annular domains, two Cauchy integrals are needed to approximate the confor- 
mal map on the interior of the annulus. One integral is defined on each boundary 
of the canonical annulus. These integrals are approximated using series, comput 
ing the coefficients with FFT’s. The Horner summation scheme is then applied to 
evaluate the truncated series at the grid points in the computational domain. 

Figure 4 shows an annular elliptic region with outer boundary of aspect ratio 2. 
The corresponding region in the computational domain is rectangular with a - 0, 
b = tt/ 2, and c = 0.42, in 2.2. Recall that the conformal map from the rectangular 
computational domain onto the elliptic region is obtained by composing the expo- 
nential map with the Wegmann map. Figure 5 shows a thin annular elliptic region 
with outer boundary of aspect ratio 3. In this case o. 0, 6 ^r/2, and c 0 

2 . 2 . 

5.0 PRACTICAL LIMITATIONS OF THE WEGMANN METHOD 

In the preceding sections, the Wegmann methods were shown to be extremely 
powerful techniques to generate the conformal coordinate transformations for simply 
and doubly connected regions. Using FFT’s and the Horner summation scheme, 
the overall Wegmann grid generation method requires 0[N log N) computations, 
where N is the number of discretization points on the boundary of the computational 
domain. (The operation bound is based on the assumption that the total number of 
grid points is less than the value of N.) This order of computations is an improvement 
over standard integral techniques which typically require 0{N 3 ) computations. The 
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standard integral approach is to compute the inverse coordinate map, and then 
approximate the desired map using an interpolative process. As an example of an 
integral approach for the inverse map, see Symm [8]. 

The Wegmann technique is also more efficient than the numerical method de- 
veloped by Fornberg, [9], [10], and |ll]. For simply connected domains the Weg- 
mann scheme is approximately seven times more efficient than the Fornberg scheme 
(Wegmann |ll]). For doubly connected domains, the Wegmann method is even 
more efficient than this, since it converges quadratically while the Fornberg method 
converges linearly. It may be the case that the Fornberg method is easier to code 
than the Wegmann method. It may also be the case that the Fornberg method is 
easier to apply to an arbitrary region with a smooth boundary, since the Wegmann 
method is sensitive to the initial guess for the inverse boundary correspondence 

function. 

Some limitations of the Wegmann method determined during the numerical 
experiments will now be discussed. A basic limitation of the Wegmann method is 
in the choice of functions used to describe the boundaries of physical regions. In 
general, Wegmann [l, 2, 3] showed that a function with Holder continuous deriva- 
tives may be used to generate the inverse boundary correspondence function needed 
to construct the desired conformal mapping. However, if simple smooth polar ex- 
pression s of the form: 

z{6) = p{8) exp(i'fl), 5 - 1 


are used to represent the boundary of the physical region, the Wegmann method 
does not necessarily converge. The difficulty is that, although the first three deriva- 
tives of z are Holder continuous, they may become large. This forces the initial 
guess for the Wegmann technique to be very good and, in general, a very good 
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initial guess is not available. 

The polar form of the boundary functions has been implemented for simply and 
doubly connected elliptic regions. Several initial guesses based on linear and non- 
linear functions for the inverse boundary correspondence function have been tried 
in conjunction with the polar representation. In all cases the Wegmann method 
fails to converge. 

In order for the Wegmann method to converge without a good initial guess, 
holomorphic (complex, analytic) functions are used to describe the boundaries of 
the physical region. In the case of the simply and doubly connected elliptic regions 
the map: 

asin (— 0 + ip) 5.2 

has been used to represent the boundaries of the regions. Here a and p are real 
numbers and the variable 6 is defined in the interval (0, 27r]. Figure 6 shows the 
first derivative with respect to 6 of the holomorphic and polar representations for an 
aspect ratio 2 ellipse. Clearly, the derivative of the holomorphic function does not 
fluctuate as much as the derivative of the polar representation. It should be noted 
that the holomorphic representation must satisfy the Cauchy-Riemann equations 
whereas the polar representation does not. 

In contrast to the polar form, the holomorphic form with an initial guess of: 

t{6) = 0, 5.3 

almost always converges. For the simply connected method, regions with aspect 
ratio up to 4 have been run successfully. 

In the case of the doubly connected method, an initial guess for the inner 
radius, /i, of the canonical annulus is also required. When a good guess for \i is 
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given the doubly connected Wegmann method works very well. However, for a case 
where the outer ellipse had an aspect ratio of 2 and the inner ellipse had an aspect 
ratio of 10 the scheme failed to converge. 

In an attempt to improve the performance of the Wegmann me thod in such 
cases an intermediate conformal mapping, J, is introduced to precondition the 
physical space. Let P denote the doubly connected physical region. Let 1 denote 
the image of P under the map J. The goal is to use the Wegmann method to 
construct the conformal map from a canonical annular region, C, onto J, and then 
apply the inverse map J~ 1 to obtain the desired conformal map from the canonical 
domain onto the physical region. This preconditioning has been used successfully. 
For the present example, the preconditioning map is the inverse Joukowski map 
defined by 



where 7 is a function of the dimensions of the inner elliptic boundary. The inverse 
Joukowski map carries the inner elliptic boundary onto a circle. The outer ellipse is 
mapped onto a curve which is closer to a circle than the initial outer curve. Figure 7 
shows the image of the annular region under the inverse Joukowski map. 

The second limitation of the Wegmann method involves the number of dis- 
cretization points needed for regions bounded by high aspect ratio curves. The 
problem of determining a conformal mapping numerically generally leads to an ill- 
posed computation. One cause of this is the local angle preservation of conformal 
maps. In simply connected regions bounded by high aspect ratio curves, the local 
orthogonality forces evenly spaced discretization points in the computational do- 
main to be crowded together in the given region. The crowding phenomenon causes 
the ill-posed numerical properties seen in computing the conformal maps. 
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Figures 8 and 9 show the crowding of the boundary discretization points for 
the simply connected elliptic regions of aspect ratio 2 and 3, respectively. Clearly, 
as the aspect ratio of the boundary curve for a simply connected domain increases, 
the crowding phenomenon becomes more severe. 

The Wegmann method works well in determining the solution of the boundary 
correspondence problem for simply connected elliptic regions. However, as the 
aspect ratio increases, the number of terms needed in the series approximation of 
the desired conformal map also increases. Table 3 shows the maximum error found 
on the boundaries of the ellipses as a function of aspect ratio, and the number of 
terms used in the series approximation. The error shown in Table 3 is defined by 
evaluating the series approximation of the conformal map at points on the boundary 
of the canonical domain and substituting the result into the expression. 

Error - |F(x,y) - l|, 

where F(x,y) is defined by, 



Table 3 shows that a large number of terms are required in the approximating series 
as the aspect ratio of the domain increases. Therefore, these error results suggest 
that the Wegmann grid generation technique is practical to apply to simply con- 
nected elliptic regions of aspect ratio less than 4. For elliptic regions of aspect ratio 
greater than 4, it may be possible to apply preconditioning to reduce the number of 
terms required in the approximating series. However, this form of preconditioning 
has not been attempted in the present study. 

The severe crowding phenomenon exhibited by the conformal maps for simply 
connected regions has not been observed for the doubly connected regions. The 
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elliptic annular regions for which the Wegmann method has been tested show little 
or no crowding effects. Figure 10 shows the image of the discretized points for 
a typical annular elliptic region. The lack of crowding suggests that the annular 
elliptic regions may be approximated with fewer terms than required by a simply 
connected region of the same outer aspect ratio. 

6.0 SUMMARY 

This paper has presented numerical experiments for the implementation of the 
mapping techniques of Wegmann. Through the use of FFT’s this method requires 
0(N log N) operations. In addition, the technique is quadratically convergent. Sev- 
eral examples for both simply and doubly connected regions have been examined. 
In particular elliptic regions have been considered. The Wegmann techniques have 
been shown to work for high aspect ratio elliptic regions. However, the mam prac- 
tical limitation of these methods for both simply and doubly connected regions has 
been found experimentally to be the functional form used to represent the bound- 
aries. For general regions, the Wegmann method works well if the boundaries are 
represented by holomorphic functions. However, the methods may not converge 
at all if a general smooth polar representation is used without a very good initial 

guess. 

In spite of these limitations the Wegmann method is more efficient than other 
conformal mapping techniques. A further application of this scheme is in the de 
velopment of conformal maps for single and double elements airfoils. For this case 
it is necessary to compute the maps on the exterior of a domain. A comparison 
between the Wegmann technique and other existing methods for this problem 
being undertaken by the authors. 
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k 

N = 64 

N = 128 

N = 256 

A 7 = 512 

1 

2.6(-2) 

2.6(-2) 

2.6(-2) 

2.6(-2) 

2 

2.7 (-4) 

2.8(-4) 

2.8(-4) 

2.8(-4) 

3 

8.0(-8) 

3.2(-8) 

3 .2 (-8) 

3.3(-8) 

4 

9.4(-9) 

6.0(-14) 

4 .5 (- 1 5) 

3.5(-15) 

5 

1.6(-9) 

1.0(-14) 

8 .5 (-1 5) 

6.5(-15) 


Table 1: Error results for the inverse boundary correspondance function of the inverted 
ellipse based on the simply connected Wegmann scheme. 


k 

N = 64 

A 7 = 128 

N = 256 

A 7 = 512 

1 

2 

3 

4 

5 

2.2(-l), 1 • 8(- 1 ) 
1.4(-2), 6.9(-3) 
8.3(-5), 4.l(-5) 
2.1 (-9) , 2.5 (-7) 
1.4 (-9), 3.4 (-7 ) 

2.2(-l), l-8(-l) 
1.4(-2), 6.9(-3) 
8.4(-5), 4 . 1 (-5) 
2.1 (-9) , 2.6(-7) 
1.4(-9), 3.6(-7) 

2.2(-l), 1.8(-1) 
1.4(-2), 6.9(-3) 
S.4(-5),4.l(-5) 
2.2(-9), 2.6(-7) 
1.4(-9),3.6(-7) 

2.2(-l) , l-S(-l) 
1.4(-2), 6.9 (-3) 
8.4 (-5) , 4.l(-5) 
2.2(-9), 2.6(*7) 
1.4(-9), 3.6(-7) 


Table 2: Error results for the boundary correspondance functions of the non-concentric 
annular region based on the doubly connected Wegmann scheme. 


N 


■sm 

AR = 4 

16 

1.68(-l) 

NC 

NC 

32 

3.52(-2) 

NC 

NC 

64 

4.49(-3) 

2 .2 1 (- 1 ) 

NC 

128 

1.62 (-4) 

1 . 1 8(- 1 ) 

NC 

256 

1.95(-7) 

4.68(-2) 

NC 

512 

1.30(-12) 

7.31 (-3) 

1 -8 1 (- 1 ) 

1024 

NC 

1.06(-3) 

i.io(-i) 

2048 

NC 

4.81 (-5) 

5.60(-2) 


Table 3: Maximum error found on the boundaries of simply connected regions of aspect 
ratio AR as a function of the number of terms in the series approximation. NC 
denotes not computed. 
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Figure Captions 


Figure 1. An example grid for a simply connected region bounded by an aspect 
ratio 2 ellipse. 

Figure 2. An example grid for a simply connected region bounded by an aspect 
ratio 3 ellipse. 

Figure 3. The computational space for the Wegmann maps for the aspect ratio 2 
and 3 ellipses. 

Figure 4. An example grid for a doubly connected region with an outer boundary 
defined by an aspect ratio 2 ellipse. 

Figure 5. An example grid for a doubly connected thin region with an outer bound- 
ary defined by an aspect ratio 3 ellipse. 

Figure 6. Comparison of the derivatives of the boundary representations. The 
holomorphic case is denoted by the dot— line curve. The polar case is denoted by 
the dotted curve. 

Figure 7. The image of an annular region with an elliptic outer boundary of aspect 
ratio 2 and an elliptic inner boundary with aspect ratio 10, under the action of the 
J map. 

Figure 8. Plot of the IBCF for the aspect ratio 2 ellipse with N 64. 

Figure 9. Plot of the IBCF for the aspect ratio 3 ellipse with N 64. 

Figure 10. Plot of the IBCF for the ellipses of aspect ratio 2 and 3 with N = 64. 
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Figure 1. An example grid for a simply connected region bounded by an aspect 
ratio 2 ellipse. 



N 

X 
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Figure 4. An example grid for a doubly connected region 
defined by an aspect ratio 2 ellipse. 
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Figure 7. The image of an annular region with an elliptic outer boundary of aspect 
ratio 2 and an elliptic inner boundary with aspect ratio 10, under the action of the 
J map. 
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ABSTEACT 

Thi* paper describe* a calculation technique for de- 
termining the stability of jet. of arbitrary cross section. 
In particular, elliptic and rectangular jet* are con.id- 
ered The numerical procedure involves both » 
formal transformation between the computational do- 
main and the physical plane and a solution of the trans- 
formed stability equation in the computational domain. 
Modern, efficient, conformal mapping* are used for bo 
simply and doubly connected domains. The numeric* 
solution is based on a finite difference/pseudospectral 
discretiration of the stability equation. The technique 
is verified by comparison with previous calculations for 
circular and elliptic jet*. Calculations are performed for 
the stability of elliptic and rectangular jets of aspect 
ratio 2. Growth rates, phase velocities, and pressure 
eigenfunctions are calculated. 


1. INTRODUCTION 

This study is motivated by the authors’ interest in 
turbulent mixing in free shear flows. It is now generally 
acknowledged that the mixing process is dominated by 
the dynamics of large scale coherent structures. In ad- 
dition, the local properties of these structures may be 
modeled by a linearised analysis. This has been demon- 
strated in the experiments and analysis of excited free 
shear layers and jets by Gaster, Kit and Wygnanski 
and Petersen and Samet 2 . Tam and Morris 3 '* made use 
of instability wave model* of the large scale structures to 
predict the noise radiation from supersonic shear layers 
and the development of excited jets. In addition, Lion 
and Morris and Giridharan 6 have developed Reynolds 
stress closure schemes in which the unknown turbulent 
stresses are described by solutions of the local stability 
equation. 

The present paper is concerned with the instabil- 
ity of jets of arbitrary cross sections. This analysis is 
an essential component in the extension of some of the 
analyses of turbulent flows described above to more com- 
plex geometries. Non-circular jets have been observed 

* NASA Graduate Student Researcher, Depart- 
ment of Aerospace Engineering, Student Member, AlAA 

. Professor, Department of Aerospace Engineering, 
Member AlAA. 


to have enhanced mixing propert.es over circular jets. 

This make, their use attractive as injectors in combus- 
tors Schadow et al 7 demonstrated this improved mixing 
using an elliptical jet in a dump combustor. Rectan- 
gular jets have applications in thrust-vectonng/thrust- 

reversing engine nossles for future fighter aircraft. 

The stability of elliptic jet. has been studied by 
Crighton* and Morris®. In the former case a vortex 
.beet representation of the jet flow was used. In the 
latter case more realistic, finite thickness, shear layers 
were considered where the analytic description of the 
mean flow enabled separable solutions of the stability 
equation to be found in elliptic cylindrical coordinates^ 
Koshigoe and Tubis 1011 used both a finite element and 
finite difference approach to consider the stability of jets 
of elliptic and triangular cross sections. 

In the present paper we describe efficient algo 
rithm. to compute the stability of non-circular jets. 
These techniques are applied to elliptic and rectangu- 
lar jets. These results are a partial summary of a doc- 
toral dissertation of one of the authors (RSB) . b 
the next section the numerical methods are esen e 
This involve* a conformal mapping technique for sim- 
ply and doubly connected domabs and a hybrid finite- 
difference / pseudospectral solution of the stability equa- 
tion. The stability characteristics of elliptic and rect- 
angular jets are then given b the form of the variation 
of the growth rates and phase velocities with instability 
wave frequency and eigenfunction, for the most unstable 
modes. 

2.0 ANALYSIS AND COMPUTATIONS 

A jet flow is considered issubg from a nozzle of 
arbitrary cross section. The governbg equation is given 
in Cartesian coordbatea, (x.y.z). The axis of the jet is 
aligned with the z direction and the axial mean velocity 
is denoted by W(x , »). The mean velocity components in 
the x and y direction, are neglected. This » the parxUel 
flow approximation of hydrodynamic sta ty. ■ t 

parallel flow assumption is not applied and the effects of 
flow divergence are considered, a problem results which 
require* a multiple-scales analysis. These effects are not 
considered in the present study. 

A lbear, elliptic, partial differential equation for 
the pressure is obtabed by takbg the divergence of the 
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momentum equation and using the equation of continu- 
ity. The resulting equation if linearised about the mean 
flow. The velocity fluctuation! are eliminated b favor of 
the prewure fluctuation ufbg the lbeariied momentum 
equation!. If the prewure fluctuation if written b the 

form. . . ■ _ / . . l 

p(x,y,z,t) = pix.y)* 1 , (2.1) 


where, 

7 = az — wt, 

and a u the axial wavenumber and u i* the instability 
wave frequency, then the equation for p may be written 

(A-a^p+^^VW-V^O (2.2) 

Eqn. (2.2) is the Rayleigh equation governing the invis- 
cid, incompressible, spatial stability of jets of arbitrary 
cross section. In order to determine the pressure p, in 
eqn. ( 2 . 2 ), boundary conditions must be added. In this 
case, p satisfies: 


and 


p — ► 0 as |x| — * -foe (2*3) 

p is finite as |x| — ► 0 (2-4) 


The computational goal in solving the stability problem 
defined by eqns. (2.2)-(2.4) is to determine the complex 
wavenumber spectrum for a given frequency, w. Several 
difficulties arise in computing the wavenumbers or eigen- 
values associated with the Rayleigh problem. Firstly, 
the wavenumber, a, appears nonlinearly in ( 2 . 2 ). This 
nonlinearity prevents the direct use of standard matrix 
eigenvalue calculation schemes. Secondly, if the stability 
analysis is used to model turbulent mixing in jet flows, 
the stability problem must be solved at a large number of 
cross sections in a given jet flow. Therefore, any scheme 
developed must be computationally very efficient. 


2.1 The Generalised Rayleigh Equation 

To consider jets of arbitrary geometry, conformal 
mappings are used to map standard computational do- 
mains onto realistic jet flow cross sections. Conformal 
mapping is analytically a very desirable technique, since 
these maps simplify the governing differential equation 
by generating a diagonal metric tensor. Therefore, con- 
formal maps are applied to simplify the Rayleigh prob- 
lem for a jet of arbitrary geometry. 

Recently, very efficient schemes have been devel- 
oped to determine the conformal maps from standard 
computational domains onto arbitrary regions in the 
plane. Wegmann 13,i4,ls proposed a scheme to compute 
the conformal maps from canonical domains onto simply 


and doubly connected regions with smooth boundaries. 
This technique has been applied to determine the coor- 
dinate maps needed for a realistic elliptic jet flow cross 
section. For polygonal regions, Trefethen 1617 has devel- 
oped an efficient software package, SCPACK, to deter- 
mine the conformal map and its inverse from the interior 
of a polygon onto the unit disc. SCPACK has been ap- 
plied to determine the conformal maps needed in the 
study of a rectangular jet flow cross section. The details 
of the conformal grid generation techniques are given by 
Baty 13 . To apply conformal maps, the Rayleigh problem 
must be recast in general conformal coordinates. Let 
Cartesian coordinates be denoted by (x 1 ,! 3 ), and let 
the computational coordinates be denoted by (y 1 , y 2 ) * 
Now, let / denote a conformal map satisfying the fol- 
lowing relations: 

= iM/(y l +*V)) ( 2 - 5 ) 

and 

= /m(/(y l +»y 3 )) (2-6) 

Applying eqm. (2.5)— (2.6), the Cauchy-Riemann equa- 
tion!, and the general tenaor form of (2.2) then^ yield* 
the generalised Rayleigh equation in terms of (y ,y )• 

vp ‘° (r7) 


where 


g = |/V+iy 2 )f 


2.2 The Hybrid Method 

This section develops a hybrid numerical scheme 
to solve the inviscid, incompressible stability problem 
defined by eqn. (2.2). Hybrid techniques are numeri- 
cal methods which combine series approximations with 
finite difference calculations. For two dimensional par- 
tial differential equations, such as the stability problem, 
hybrid methods generate discretisation matrices whose 
order varies linearly with the approximating series sum- 
mation bound, N. This implies that the number of op- 
erations required to compute an eigenvalue of equation 
( 2 . 2 ) is on the order of: 

0(N 3 ) (2-8) 


The estimate given by (2.8) is a great improvement 
over the number of operations required by spectral and 
pseudospectral methods. These series methods gener- 
ate discretisation matrices from two dimensional series 
approximations and produce matrices whose sise is pro- 
portional to: 

(N + l ) 3 (29) 
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Equation (2-9) implies that the* method, require a 
number of operation* on the order of: 

0((tf + l)*) (210) 


Here T (y 1 ) » the deriv * tiv * of the *** ° rder Cbeby - 
»hev ' poTynomial Evaluating (2.13) at the gnd potnt. 

leads to the relation. 


to compute an eigenvalue of the Rayle.gh problem. 
Furthermore, hybrid technique, have the advantage 
of increased accuracy over purely finite <M««nce a I> 
proache.. The accuracy of a hybnd scheme depend, on 
the propertie. of the approximating senes and on the 
accuracy of the finite difference .cheme. The pre.en 
study uses a p.eudo.pectral »erie« to develop a hybnd 
scheme. 


A p.eudo.pectral #erie. assume, that the function 
to be approximated is known, or may be computed, on 
a fixed set of point, in the computational domain. This 
information is then used with an appropriate set of basis 
functions to form a finite series approximating the func- 
tion. Basing the function’, approximation on a known 
set of points defines a grid in the computational domain. 


The hybrid method presented here is defined using 
a pseudospectral series based on the Chebyshev poly- 
nomials. The details of the pseudospectral technique 
come from the one dimensional pseudospectral theory 
in Gottlieb et al 1 ®. 


To outline the method, consider the Rayleigh 
eqn. (2.8). Assume that the function to be approxi- 
mated may be represented in a series of the form: 


( 211 > 

.=0 

Here y 1 represents the aximuthal variation, while y 2 rep- 
resents the radial variation in an arbitrary jet cross sec- 
tion. In eqn. (2.11), the coefficients <*(y,-,y 2 ). *r« func- 
tions of the grid points defined by: 


where £y, is the Kronecker delta. 

Next, substituting eqn. (11) into the Rayleigh equa- 
tion t yields: 


y>(y‘,y 3 )/;v) + iy^.y 3 )/.^ 1 ] 

* v • n 


t=0 


^ ^ ISO * 0 

-a 2 ^o(y. , ,y 3 )/.(y l )=o ( 215 ) 


» =0 


where 2a 

P = Z - aW 

To simplify this equation, the derivatives of the approx- 
imating basis functions at the grid points must be de- 
termined. Reference 18 gives these derivatives as: 

^\A(yM| = (J D p )ji (2-16) 

<%‘) p I,.-,; 

where 


c A~ 

iV+j 

,x if 

(2.17) 

" c *(yy 

-y‘) 


[D l )a = 

-yy 

2 ( 1 -(yj) 2 ) 

(2.18) 

(D 1 ) oo 

2 N* + 1 
6 

(2.19) 

(D l ) NN 

= -(I? 1 ) 00 

(2.20) 


= cos(^) for j = 0,1,2,...,* (2-12) 

and the radial direction, y 3 . Moreover, the basis func- 
tions, /*, are rational functions defined by: 

_ qi-^nwH-ir') (2 . 1S| 

" i V ■ -yj) 


with 


and 


Co = C S 


= 2 


= 1 otherwise. 


and where . * 

(£)»>) = (D l ) p (2-21) 

Now, combining eqns. (2.16)-(2.2l), eqn. (2.15) and l ap- 
plying eqn. (2.14) produces a system of linear, second or- 
dinary differential equations in terms of the coefficients: 

, 2a \ dW n j 2^ 

‘."(vi.v’l + 

- ^.(vi.v 3 ) + y’HD 3 )* 


s 


+ ( 


20 }£21£,,( S ,,'.V J )(D‘)„=0 (2.22) 


u ,-aW‘ ay 1 fz 0 
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'he eqnft. (2.22) are then recast as a collection of first 
rder systems. This collection of systems of first order 
q nations can then be recast as a matrix equation with 
he unknown vector being the coefficients of the approxi- 
lating series. The resulting matrix equation may be in- 
.ted once the boundary conditions are determined 
n added to the matrix equation. 

5 The Boundary Conditions 

In order to evaluate numerically the discretiia- 
on matrix associated with the Rayleigh equation, the 
oundary conditions must be included. In the region 
utside of the mixing layer, the velocity is constant, re- 
-icing the Rayleigh equation into the Helmholts equa* 
on: 

(A-$a 3 )p = 0 (2.23) 

Silowing Batchelor and Gill 19 , the general solution of 
-in. (2.23) in polar coordinates for the exterior region 
' given by: 

22 (A n /f^(Vor) + 5„^(»ar)) exp(Vnfl) (2.24) 

n = 0 

‘here if 1 and H 3 denote Hankel functions of the first 
od second kind respectively. Using the condition, (2.3), 
lat the pressure must approach sero as r approaches in- 
rity, implies that the boundary condition on the outer 
fge of the jet flow cross section is of the form: 

p = 21 A "^n(* Q[ r)exp(in«) (2.25) 

n = 0 

^oreover, since the pressure must be bounded as r ap- 
'*oaches sero, condition (2.4), the boundary condition 
the inner edge of the jet flow cross section becomes: 

CO 

p = ^B„J n (»ar)ex p(tn5) (2.26) 

n=0 

' here J n is the Bessel function of the first kind. Next, 
physically realisable aiimuthal terms are determined 
7 order to simplify eqns. (2.25)-(2.26). All the jet ge- 
fnetries to be computed in this study are symmetrical 
J»out both the horisontal and vertical axes in the plane. 
* !ius, from Morris 9 , the possible pressure variations in 
azimuthal direction correspond respectively to four 
^ -vases of functions depending on symmetries about the 
\jor and minor axes. This then gives the general solu- 
1 n for the physically possible boundary' solutions. The 
“ulting infinite series defining the pressure boundary 
uditions on the edges of the shear layer cross section 
come: 

oo 

22 ^nC„(tar) cos(2n5) (2.27) 

n = 0 


22 A„C„(tar) sin((2n + 1)5) (2.28) 

nsO 

f;X„C„(tar)»m((2n + 2)5) (2.29) 

n>0 

£A„C n (.or)co.((2n+ 1)5) (2.30) 

n>0 

where C n represents either J n or H\. Now, recalling 
that the jet flow cross section is assumed to be symmet- 
rical about both axes, allows the computation to be re- 
stricted to the first quadrant in the physical plane. The 
standard computational domain for this physical region 
will be a rectangle. On the edges of the computational 
rectangle which correspond to a constant radial value, 
the functions defined by eqns. (2*27 ) ( 2.30) will be ap- 
plied. Before these boundary conditions are evaluated 
in computational space, they are transformed in terms 
of the computational coordinate system using the met- 
ric generated by the conformal mapping. On the verti- 
cal edges of the computational domain, which represent 
lines of constant angle, the boundary conditions are de- 
termined by the sysmmetry conditions. If the pressure 
fluctuation is not symmetrical about an axis, that is, if 
it changes sign across an axis, then the corresponding 
boundary condition becomes: 

p = 0 (2.31) 

However, if the pressure fluctuation is symmetrical 
about an axis, that is, the sign does not change, the 
pressure boundary condition becomes: 

= o (2.32) 

36 

Notice that the boundary conditions defined by 
eqns. (2.27)— (2.30) for the horisontal edges of the com- 
putational domain are consistent with the boundary 
conditions imposed on the vertical edges of computa- 
tional domain. 

The matrix equation may be integrated explicitly in 
the radial direction once the boundary conditions have 
been converted into the appropriate initial conditions. 
The boundary conditions on the horisontal edges of the 
computational domain are converted into initial condi- 
tions using a generalised shooting method. Let A de- 
note the summation bound for the approximating series. 
Then there are N - 1 interior grid points. On the lower 
edge of the computational domain, the first term and its 
derivative from the exact series solution as given above 
are evaluated. This becomes the initial condition on the 
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lower boundary- The matrix containing the syst err. of 
differential equation, is then integrated to the geomet- 
center of the computational domain, yielding V,. 
At each »tep in the explicit integration procedure the 
boundary condition, on the vert.ca 1 ed f e. of the com- 
putational domain are satisfied by .olving for the first 
and last coefficient, of the approximating wne. or it. 
derivative. 


On the upper edge, the fir»t term in the exact ee- 
rie, solution and it. derivative are evaluated also. Then 
these value, are used to integrate the matrix equation 
to the center of the computational domain, producing 
V". This process i* repeated for each term in the series 
using N- 1 term, from the exact series solution, on the 
horiiontal edges of the computational domain. 


The resulting integrated solution, and their deriva- 
tives are then matched at the center of the domain. The 
matching is accomplished by requiring that a linear com- 
bination of the 2 (N - 1) solutions be equaled to *ero: 


£ (R.V? + S ; V') = 0 (2.33) 

»= 1 

Requiring that eqn. (2.33) have a non-trivial solution 
then forces the determinant of the matrix of integrated 
solutions to be »ero. Recalling that the solution vec- 
tors are implicit function, of a fixed real frequency and 
some guessed complex wavenumber implies that a local 
scheme may be used to determine the wavenumber*. In 
this study the Newton-Raphson scheme was used. 

Once a wavenumber or eigenvalue of the Rayleigh 
problem has been computed, its corresponding eigen- 
function may be determined. The hybrid method com- 
putes an eigenfunction as it integrates the initial condi- 
tions from the edges of the shear layer to its geometric 
center. However, the relative weights of the integrated 
solutions are not known. These weight, are precisely the 
coefficients which force the solutions and their deriva- 
tives to match in the shear layer. 

Recall that this matching is given mathematically 
by eqn. (2.33), where R< and S< are the unknown coeffi- 
cients. These coefficients may be determined by recast- 
ing (2.33) in the form: 



After the matching coefficients are determined they 
are used to scale the initial conditions. The Rayleigh 
problem is then integrated a final time using the scaled 


initial condition, and the p^udospectral ^pbtudes, 
0 (yi y»), are stored along each radial grid line There- 
fore'^ eigenfunction is approximated discretely in the 
radial direction, and by a series in the a.imuthal direc- 




2 A The Mean Velocity Profile 

The mean velocity profile used in the computa- 
tion. is based on a generalisation of the profile given by 
Michalke 20 . For the round jet Michalke chose a velocity 
profile in the mixing layer of the form: 



(2.35) 


S 

for R - ~ < r < oo 

where R is the jet radius, 6 is the momentum thickness 
and 6 is a fixed real number satisfying: 


u»i4) » i <*•») 


Notice that the velocity profile defined by eqn. (2.35) 
is only a function of the radial direction. Furthermore, 
since conformal map. are being applied to generate the 
computational coordinates, (y l »y 3 )i tiie radial coor 1 - 
nate, y 2 , will be uncoupled from the asimuthal coord)- 
nate' y», allowing a general velocity profile to b« ex- 
pressed in terms of y 2 . To generalise eqn. (2.35), let / 
denote a conformal map carrying a standard computa- 
tional rectangle onto the first quadrant of a jet fiow cross 
section. Also, let /(y 1 + »y 3 ) denote the minor axis of 
the jet cross section in terms of the computational co- 
ordinate system. Then in terms of / the generalised 
velocity profile is defined as: 


, if . J h /(v 1 ±i£l n 

"V) = 2l 1 + tanh W /(v'+^r 

for y 2 < y 2 < <«> 


(2.37) 


where y 2 is the half velocity point, B is the length of 
the minor axis, 6 B is the momentum thickness on the 
minor axis, and y 2 is a real value satisfying. 


B , 

tanh( — |i 


/( y* + t y?) n 

/( y l +»yl) 


w 


1 


(2.38) 


The momentum thickness used in the generalised mean 
velocity profile is defined on the minor ajcis by. 



-W)dy 


for 


I 


0 


(2.39) 
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Kg an example of a non-circular mean velocity profile, 
f the function aco*(y l + ly 3 ) i* used to generate the 
>hysical coordinate*, the mean velocity profile defined 
-»y eqn. (2.37) reduce* to. 



for y J < y 3 < oo. 

"he next section will outline teat results for the stability 
ode using the velocity profiles defined in this section. 

5.0 NUMERICAL RESULTS 

Ln this section the generalised Rayleigh problem 
'overning the Linear inviscid stability of incompressible 
ets of arbitrary geometry is solved using the hybrid 
cheme described in the previous section. To consider 
^ts of arbitrary geometry, conformal transformations 
<ave been used to map standard computational domains 
*nto jet flow cross sections in physical space. Calcula- 
ions are performed for the stability of the annular shear 
ayer region of rectangular and non-confocal elliptic jets 
-,f aspect ratio 2. These calculations are performed for 
»zimuthal normal modes corresponding to the flapping 
»nd varicose instabilities observed in non-circular jets, 
irt addition, examples of the eigenfunctions for these 
ases are shown. Firstly, the validation of the stability 
-ode, by comparison of its results to benchmark calcu- 
lations for the round and confocal elliptic geometries, is 
riven. 

t.l Code Verification 

The stability code has been validated numerically 
Sr several different geometries and boundary condi- 
J *ons. The series of numerical tests performed involved 
imputing the eigenvalues associated with the maxi- 
mum rate of growth for the flapping and varicose ai- 
''nuthal normal modes of incompressible circular and 
‘ trifocal elliptic jets. These results have been compared 
the results of Morris 9 , who computed the wavenum- 
bers associated with the maximum growth rate for in- 
impressible confocal elliptic jets. 

In order to compare the present computation with 
Morris 19 results, a relative error is introduced. Define: 

% Error = (3.1) 

\*b\ 

' here a p are the current results and a are the published 
"suits. 

For the circular case, the complex exponential map 
'* used to generate the grid for the generalized Rayleigh 


problem. The mean velocity profile used in this com- 
putation is given by eqn. (2.35) above. Ln all the test 
case* the momentum thickness on the major axis, B At 
is taken to be 0.02, as in ref. 9. Also, all the computa- 
tions assumes geometries such that the velocity profile 
satisfies: 

.001 < VV < .999 (3.2) 

in the shear layer. Table 1 compares the wavenumbers 
computed using the hybrid technique, for the aspect ra- 
tio 1 case, with those of ref. 9 for the aspect ratio of 
1 . 001 . 


( Mode 

1 Fr*q 

Hybrid ■ 

Morru 

Error j 

i V 

5.441 

1 10.20T-5.6S7 i j 

10. 199-5. 6S5i 

0.07 

1 f a 

5.453 

! 10. 244-5.6521 

10. 238-5.6491 

0.06 1 

Fb 

5.456 

10.253-5 652i | 

lC.243-5.651t 

i 0.08 | 


Table 1 Aspect ratio 1.0 results for the hybrid 
method compared to previous calculations 9 . 

In Table 1, F A and F B represent aiimuthal flapping 
modes about the major and minor axes respectively, 
while V represents the varicose aiimuthal mode. The 
results of the stability code are approximately indepen- 
dent of the number of collocation points in the aiimuthal 
direction. For the circular case, a Runge-Kutta integra- 
tion scheme is used with a fixed step-siie of 0.004. 

For the confocal elliptic case, the complex cosine 
function is used to calculate the gnd for the general- 
ised Rayleigh equation. This coordinate transformation 
generates elliptic cylindrical coordinates in the physical 
plane, as used in ref. 9. The confocal elliptic thin shear 
layer calculations have been performed for the aspect 
ratio 2 case. The mean velocity profile used in both cal- 
culations is a special case of the general profile of the 
last section and is given by eqn. (2.40). Table 2 com- 
pares the results of the hybrid method to the results of 
ref. 9 for the aspect ratio 2 confocal elliptic shear layer. 


Mode 

Fre<; 

Hybrid 

Morn* 

y^Erro: 1 

V 

I 5.657 

10. 156-4.496* 

10.1 35-4.507i 

0.21 

f a 

5.010 

9.307-3.6641 

9.322-3.677i 

0.19 

Fb 

j 5.657 

10.045-4.472i 

10.027-4.507i 

0.35 


Table 2 Aspect ratio 2 results for the hybrid method 
compared to previous calculations 9 . 

The step size is 0.006 for the aspect ratio 2 case. For 
the aspect ratio 2 confocal elliptic shear layer, the results 
shown are computed with 7 interior collocation points. 
If less than 4 interior collocation points are used the 
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computed wavenumber, exhibit a large error in com par- 
uon with the benchmark value*. Conversely, if a large 
number of interior collocation point, are used, .ay above 
9 it become, very difficult to locate the wavenumber, 
associated with the maximum growth rate. By adding 
more interior point, to the computation, two distinci 
limiting ^procewes are affected. Firstly, the approxima- 
tion of the eigenfunction become, more accurate because 
more geometric information about the physical domain 
b supplied to the approximation. Secondly the func- 
tions defining the boundary conditions on the edges o 
the shear layer become more accurate by including erms 
which fluctuate more rapidly in the aiimuthal direction. 
It is believed that adding solutions which represent t e 
rapidly fluctuating aiimuthal terms to the integrated 
discretisation leads to a determinant minimuation prob- 
lem which is ill-conditioned. Presently, this represents 
the main difficulty found in solving the Rayleigh pro 
lem with the hybrid method. Since the generalise 
shooting technique couples the geometry and the bound- 
ary conditions, one possible way of correcting the prob- 
lem would be to normalise appropriately the determi- 
nant equation before it is minimised. No further work 
on this problem has been attempted in the present study. 


In the analysis of the instability of jets of arbitrary 
geometry the eigenfunctions associated with the com- 
puted wavenumbers may also be determined. Recall 
that in this case the eigenfunctions correspond to the 
pressure function, p. All the other field variables may 
be related to this function. Thus, once the pressure 
eigenfunction is known, the distributions of the velocity 
components associated with instability waves or large 
scale coherent structures are determined completely. 


The technique for the evaluation of the eigenfunc- 
tions was described in Section 2.3. Two numerical 
checks have been performed in analy.ing the results of 
the eigenfunction calculation for the round jet. Firstly, 
the stored integrated solutions have been shown to 
match at the geometric center of the computational 
domain. Secondly, the computed constants, given in 
eqn. (2.34) which weight the initial conditions have been 
checked to verify that the code is predicting the funda- 
mental instabilities for the varicose and flapping modes. 
For both of these numerical tests the code performed 
well. Further numerical checks have shown that the code 
also predicts correctly the higher order modes. 

As a verification of the eigenfunction calculation 
procedure the eigenfunctions for the varicose (axisym- 
metric) and flapping (helical) instabilities of the round 
jet have been computed. The corresponding eigenvalues 
are given in Table 3. The eigenfunctions are shown in 
Figs. 1 and 2 in the form of iso-pressure contours. The 


contour, are plotted for the fundamental varicose and 
flapping instabilities given by eqn. (2.1) with equal 
to 0. All the contour, shown are for positive value, as 
there are no negative value, for the aero phase case^ Be- 
low, when the eigenfunction, are shown for the elliptic 
and rectangular ca~s the positive and negative contours 
are shown separately. These plot, further assume hat 
the pressure field is normalised by the modulus of us 


Mode 

FVeq 

Wavenumber | 

V 

Fa 

5.44 

5.45 

10. 20-5.681 ~1 
| 10. 24-5.65i | 


Table 3 The frequencies and wavenumbers used to 
compute the eigenfunction, associated with the round 


The eigenfunction for the fundamental vancose 
mode is shown in Fig. 1. lu this case the derivative of 
the pressure with respect to the asimuthal direction is 
sero on both axes. Moreover, as expected, the pre«ure 
contours show little aiimuthal variation. The minor as- 
imuthal variation in the pressure eigenfunction indicates 
that the eigenfunctions contain very small contributions 

r 1 •_ 1 



Fig 1 lso-pressure contours for the varicose mode: 
circular jet case, T = 0. (Positive contours only shown) 

The contour plot for the fundamental flapping in- 
stability about the horisontal axis is shown in Fig. 2. For 
this mode, the pressure contours should be symmetries 
about the vertical axis and »ero along the horuonta 
axis. For both of these examples, the hybrid method 
works well in predicting the properties of the normal- 
ised pressure shape function. 
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? ig. 2 Iso-pressure contour* for the flapping mode: 
-irculax jet case, 7 = 0. (Positive contour* only shown) 

3.2 Stability Calculations 

For the numerical experiments, a non-confocal el- 
iptic shear layer and a rectangular shear layer of aspect 
*atio 2 arc chosen. Both of the jet shear layers are non** 
*1 imens ion ali* ed by requiring that the semi-major and 
^emi— minor axes, A and B, satisfy: 

| = 2 (3.3) 

md 

\TaB = 1 (3.4) 

\long the half velocity line in the shear layer. Further- 
'nore, both cases assume the hyperbolic tangent mean 
'elocity profile defined by eqn. (2.37). All numerical 
‘ests have been performed assuming a constant mo- 
mentum thickness of 0.02 in the computational domain. 

the physical plane this results in a non— uniform ax- 
'rnuthal variation in the momentum thickness: the mo- 
mentum thickness on the minor axis being greater than 
'hat on the major axis. This corresponds to the experi- 
mental profiles observed by Seiner et al 21 for an elliptic, 
Aspect ratio 2, supersonic jet. 

The first step in analysing the instability of the 
•ion— circular jets is to compute the complex wavenum- 
bers or eigenvalues for a set of real frequencies. The 
1 ases considered are the varicose and flapping aximuthal 
‘nodes. The flapping case is for flapping about the 
‘Major axis. The wavenumber* are computed by fixing 
» frequency and then making an initial guess for the 
Wavenumber. A Newton-Raphson iterative scheme is 
v »sed to locate the eigenvalues. The computed wavenum- 
bers and their corresponding frequencies are used to de- 
’ ermine the local growth rate and phase velocity for the 
^stability waves. 


For the shear layer* considered, the hybrid method 
was run with both 5 and 7 interior collocation points. 
The difference between wavenumber* for these cases u 
typically in the second and third decimal places. All 
the calculations for the elliptic shear layer are based on 
7 interior collocation points, while those for the rect- 
angular shear layer are based on 5 interior collocation 
points. Figure 3 shows the variation of the axial growth 
rate with frequency for the varicose mode of the aspect 
ratio 2 elliptic jet. The maximum growth rate is slightly 
lower than that determined for the confocal elliptic shear 
layer: see Table 2. The variation of the phase velocity, 
given by w/a r , for this case is shown in Fig. 4. This re- 
sult is typical of all the calculation* for both the varicose 
and flapping instabilities in the elliptic and rectangular 
jet cases. For the varicose instability there is generally a 
slight decrease in the phase velocity at low frequencies. 
However, in all the cases considered, the phase velocities 
of the instability waves are approximately 60 percent of 
the centerline velocity. These results are in agreement 
with those of Koshigoe and Tubis 10,11 . 


4.0 n 



2.0 4 . 0 6.0 6-0 

U) 


Fig. 3 Variation of the axial growth rate with fre- 
quency. Elliptic jet, aspect ratio 2, varicose mode. 

To determine the most unstable mode the three 
largest growth rates were interpolated using a second 
order polynomial. Table 4 shows the frequencies of the 
maximum growth rates for the elliptic jet. 


Mode 

FYeq 

Wavenumber 

V 

Fa 

4.49 

4.16 

7.86-3.5CH 

7.75-2.95i 


Table 4 The frequencies and wavenumbers used to 
compute the eigenfunctions associated with the non- 
confocal elliptic jet. 
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Fig. 4 Variation of the phase velocity with frequency. 
Elliptic jet, aspect ratio 2, varicose mode. 

Figure 5 ahows the variation of the axial growth 
rate with frequency for the varicose mode of the aspect 
ratio 2 rectangular jet. The maximum growth rate is 
much lower than that for the elliptic jet. In addition, 
the frequency for the maximum growth rate is also re- 
duced. Table 5 shows the frequencies for the maximum 
growth rates for the varicose and flapping instabilities. 
This frequency should give an indication of the initial 
vortex shedding frequency for the jet. It is not clear 
whether the calculated reduction in this frequency for 
the rectangular jet is due to the change in the geometry 
or to the distribution of momentum thickness around 
the jet. This question is being addressed by the authors 
in additional calculations. 



0.0 2.0 4.0 6.0 B.O 10.0 


Mode 

Freq 

Wavenumber 

V 

ft J 

3.16 

2.90 

5.71-1.92i | 

5.46-1.98t J 


Table 5 The frequencies and wavenumbers used to 
compute the eigenfunctions associated with the rectan- 
gular jet. 

ure 6 shows the iso-pressure contours for the most un- 
stable varicose instability in the elliptic jet and Fig. 7 
shows the corresponding contours in the rectangular jet 
case. Positive contours are shown in the upper half of 
the figure and negative contours are shown in the lower 
half. The sero contour appears in both. The most 
notable feature in both figures is the lack of regular- 
ity, compared to the circular jet contours shown above. 
The pressure fluctuations are confined to several regions 
with no apparent relationship to the particular geome- 
try. Very similar distributions are found for the flapping 
modes so that they are not shown here. 



Fig. 5 Variation of the axial growth rate with fre- g ko-pressure contours for the elliptic non- 

quency. Rectangular jet, aspect ratio 2, vancose mode. coaf(Xil j et Aspect ratio 2, varicose mode, y = 0. (a) 

positive contours, (b) negative contour*. 

As a final calculation we consider tbe pressure 
eigenfunctions for the elliptic and rectangular jets. Fig- 
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Fig. 7 Iso-pressure contours for the rectangular jet. 

Aspect ratio 2, varicose mode, 7 = 0. (a) positive con- 
tours, (b) negative contours. 

Though the contours of equal pressure level show 
no apparent structure this is not the case for the con- 
tours of equal Reynolds stresses and their gradients. In 
the rectangular jet case the distributions indicate that 
the jet is developing initially as two independent two- 
dimensional shear layers. Whereas, in the elliptic jet 
case there is a continuous variation from the major to 
the minor axes with the dominant fluctuations, for the 
modes considered being close to the major axis. Fur- 
ther details axe contained in ref. 12 and will be reported 
elsewhere. 

4.0 Summary 

This paper has presented a calculation procedure to 
determine the stability of jets of arbitrary cross sections. 
Calculations have been performed for elliptic and rect- 
angular jets. Both eigenvalues and eigenfunctions for 
the pressure have been calculated. The velocity com- 
ponents may be obtained from the pressure using the 
linearized equations of continuity and momentum. In 
turn this enables the second-order statistics, including 
the normal and shear stresses, associated with the in- 


stability waves to be calculated. If it is argued that the 
mixing process in free shear layers is dominated by large 
scale structures and that, locally, they may be modeled 
as instability waves, these second-order statistics are all 
that is needed to provide a turbulence closure scheme. 
This technique is presently being developed by the au- 
thors. 

The calculations presented here are for incompress- 
ible flow. They are readily extended to the compressible 
flow case and examples of this calculation are contained 
in ref. 12. in addition, this reference contains details 
of the conformal mapping techniques as weLl as calcu- 
lations of the velocity, Reynolds stress, and Reynolds 
stress gradient distributions. These calculations are be- 
ing used by the authors to help understand the axial 
development of non— circular jets. The results of this 
analysis will be presented later. 

This work was supported by NASA Langley Re- 
search Center under NASA Grant NAG-1-357. The 
technical monitor is Dr. J. M. Seiner 
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ABSTRACT 

The shock structure in non-circular supersonic jets 
is predicted using a linear model. This mode! includes 
the effects of the finite thickness of the mixing layer and 
the turbulence in the jet shear layer. A numerical solu- 
tion is obtained using a conformal mapping grid genera- 
tion scheme with a hybrid pseudo— spectral discretisation 
method. The uniform pressure perturbation at the jet 
exit is approximated by a Fourier-Mathieu series. The 
pressure at downstream locations is obtained from an 
eigenfunction expansion that is matched to the pressure 
perturbation at the jet exit. Results are presented for 
a circular jet and for an elliptic jet of aspect ratio 2.0. 
Comparisons are made with experimental data. 

1. Introduction 

Broadband shock associated noise is one of the ma- 
jor components of the noise of supersonic jets operating 
at off-design conditions. As a result, in recent years, ef- 
fort has been focussed on understanding the characteris- 
tics and the generation mechanisms of broadband shock 
associated noise. The early work of Harper-Bourne and 
Fisher 1 has been followed by a number of experimental 
and theoretical studies on this topic. Investigations have 
been conducted by Tanna. 2 , Seiner and Norum 2, *, Seiner 
and Yu 5 , Norum and Seiner 6, ' and Tam and Tanna 8 
among others. The noise generation mechanism pro- 
Dosed by Tam and Tanna 8 was used by Tam 9 to develop 
a stochastic model for broadband shock associated noise. 
He obtained very good agreement with experiments for 
.he near and far field noise spectra and directivity. In 
a is study the shock cell structure was modelled using the 
method of multiple-scales of Tam, Jackson and Seiner 10 . 

Recently there has been considerable interest in 
ion-circular supersonic jets with a view of achieving bet- 
er mixing characteristics and a reduction in radiated 
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noise. There has been some progress in the develop- 
ment of supersonic jet noise theories for non- circular 
jets. Tam 11 used a vortex-sheet model for the jet and 
predicted the screech tone frequencies in rectangular and 
elliptic jets. Morris, Bhat and Chen 12 used a bound- 
ary element method to predict the shock cell structure 
and screech tone frequencies in jets of arbitrary geome- 
try. Once again, a vortex-sheet model was used to de- 
scribe the jet. Morris and Bhat 13 extended their anal- 
ysis of non-circular geometry jets to include the effects 
of finite mixing layer thickness using realistic and con- 
tinuous mean velocity and density profiles. They also 
included the dissipative effects of the small-scale turbu- 
lence through the addition of eddy-viscosity terms in the 
momentum equations. However, they encountered some 
convergence problems with their numerical technique. 

In this paper, the shock cell structure of non- 
circular supersonic jets is modelled using a linearized 
analysis. In the present study a new numerical scheme is 
introduced. This method uses conformal mapping with 
a pseudo-spectral hybrid discretization scheme. This 
work is an attempt by the authors to develop models 
with a robust numerical scheme to predict the shock 
cell structure in an efficient manner. The physical model 
used here is similar to the one used earlier by Tam, Jack- 
son and Seiner 10 and Morris and Bhat 13 . In the latter 
case a body-fitted coordinate system was used to set up 
the problem for non-circular geometry jets. However, 
this scheme had convergence problems. The new nu- 
merical method uses the conformal mapping technique 
developed by Wegmann 14 and used by Baty 15 in his 
analysis of the inviscid instability of arbitrary geometry 
jets. This method transforms the physical space to a 
rectangular computational domain using a series of con- 
formal mappings. In the computational space, the flow 
variables are represented in a series form using a hybrid 
pseudo-spectral approximation. 

The next section outlines the development of the 
governing equations in terms of a general conformal co- 
ordinate system. A description is also given of the im- 
plementation of the conformal mapping technique. The 
discretization of the model differential equations using a 
hybrid pseudo-spectral scheme is also described. Section 
3 describes the application of these numerical techniques 
to the shock cell problem. In section 4, predictions are 




presented for a circular jet and an elliptic jet of aspect 
ratio 2.0. Also, comparisons are made with experimen- 
tal data. 

2. ANALYSIS 

The cross-section of the initial region of a jet con- 
sists of three regions: the potential core of the jet, the 
annular mixing region and the ambient fiuid surround- 
ing the jet. In the potential core and in the ambient 
fluid, the solutions to the Unearned governing equations 
can be obtained analytically. However, in the annular 
mixing region, where the mean velocity and density of 
the jet vary, the solution must be obtained numerically. 

The model used here for the shock cell structure is 
that developed by Tam, Jackson and Seiner 10 for circu- 
lar jets. This model takes into account the finite thick- 
ness of the mixing region and the effects of turbulence in 
the jet shear layer. A finite-difference technique using a 
pseudo-spectral hybrid scheme in conjunction with con- 
formal mapping is used here to extend the model to a 
study of the shock cell structure in non-circular super- 
sonic jets. 

2.1 Governing Equations 

The non-dimensional linearized governing equa- 
tions for the shock cell structure, in Cartesian tensor 
notation, are given by: 
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The reference scales are ry, the equivalent radius of the 
fully expanded jet and £/y, py, and py, the fully expanded 
jet velocity, density and pressure, respectively, p and W 
are the mean density and axial velocity of the jet. The 
quantities with a superscript S correspond to the per- 
turbations associated with the shock cell structure. A/y 
is the fully expanded jet Mach number. The assump- 
tions made in deriving the linearized equations given 
above are discussed in detail by Bhat 16 . The mean ve- 
locity is assumed to be known, either from experiments 
or predictions. It is also assumed that the mean flow is 
independent, locally, of the axial distance. The turbu- 
lent Reynolds stresses in the momentum equations have 
been modelled using a simple eddy viscosity model. The 
turbulent Reynolds number is given by, Re = UjTj/v t 
where is the turbulent eddy viscosity. The shock cell 
structure is modelled as spatially periodic waves that 


are time-independent. Using the locally parallel flow 
approximation, the perturbation quantities associated 
with the shock cell structure can be represented in a 
separable form with a periodic variation in the axial di- 
rection. This is given by, 


/ 5 (x i, 12 , 13 ) = /(ii i X2)exp(tQX 3 ) (2.4) 


where f s is any flow perturbation quantity and a is the 
complex axial wavenumber. 

The linearized eqns. (2.l)-(2.3) are written in terms 
of a Cartesian coordinate system. These equations must 
be transformed in terms of a general, nonsingular curvi- 
linear coordinate system, (y,,y 2 ) in the cross-sectional 
plane. This coordinate change can be performed by ex- 
pressing the differential operators in (2.1)-(2.3) in gen 
era! tensor form. These are obtained using the ap- 
proach of Eiseman 17 , with the axial coordinate given 
by, x 3 = y 3 . The shock cell structure equations in a 
general coordinate system are then given by, 
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where, g = det(g.y) and g.y is the metric tensor. The 
metric tensor and its inverse g' 1 are defined as, 
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and 


g i} = l/?.y ( 2 - 10 ) 

The transformation of eqns. (2.5)-(2.7) in physical 
space to a simple computational domain is accomplished 
using a inverse conformal mapping technique. The lor- 
mulation is similar to that given by Baty 16 and is out- 
lined briefly here. Consider a computational domain R, 
in the complex plane (yi.ya)- Let / be an analytic map- 
ping which maps S onto a given jet flow cross-section 
in the physical space (x,,x 2 ). It is assumed that / has 
a nonzero derivative, i.e. 

/'(*) = /'(yj +’V2) r 4 0 for a11 zeR * 2 ' n ) 
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The relation between the coordinates ( y i , 
i^xj) generated by f(z ) is given by, 

, ya) 

/(*) = }{v i + *&) = 21 + tl2 

(2.12) 

7 he analytic mapping / satisfies the Cauchy-Riemann 
quations, 

ax, dx 2 _ ^2 

3y, ~ dy 2 dy 2 ay, 

.'he derivative of / can then be written as, 

( 2 . 13 ) 
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a, rj and $ are finite real number* and 0 ia positive. 

The discretisation for the finite-difference scheme, 
discussed below, involves an approximating function de- 
fined on the interval | — 1, l) in the*aiimuthar direction, 
Vl . This requires a linear coordinate transformation of 
the form 

y, = Ai l -+B, and yj = *2 ( 2 - 21 ) 

with the transformed rectangular computational domain 
defined by: 

k = |(« lf » a ) 6 It 5 : -1 < si < 1,0 < *2 < P) ( 222 ) 


he components of the metric tensor, g x} are given by, 
Pn = 022 " \!'[ z )?> 9 n = 02i = 0 ( 2 . 15 ) 

7he components of g ' 3 , the inverse of g x)) are given by, 
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Then, the resulting equations for the shock cell 
iructure in general conformal coordinates can be writ- 
en, 
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nere A is the standard Cartesian Laplacian in terms 
the coordinates (yi,^)* The conformal mapping is 
rnerated numerically and the method is described in 
.e next section. 


.2 Conformal Mapping 

Let P denote a given jet flow cross-section in the 
rtysical space and C be the unit annulus with some in- 
nor radius, /x, such that 0 < p <1. Let W be the 
nformal mapping such that C is mapped onto P . As 
e desired computational domain is rectangular, an- 
ner mapping must be used to map the rectangular 
2ion Z onto the annular region C. This is given by the 
.oonential mapping exp(ir). The rectangular compu- 
'.ional domain in coordinates (yi , y?) is given by, 

< = 1(01)02) 6 R 2 : a < yi < r?,0 < y 2 < £j- (2.20) 


Computational Domain, 

t 

i 

Linear Map 

i 

Rectangular Domain, 

71 

exp(ir) Map 

! 

1 

Annular Region, C 

Wegmann Map. t 

i 

W 

Jet Cross-Section 


Physical Domain, 

V 


Fig. 1 Schematic of the transformation from the 
computational domain £ to the physical space P 

A and B may be determined once the dimensions of 
the rectangular domain Z axe known. The various steps 
involved in mapping the rectangular domain & onto the 
physical space P are shown schematically in Fig. 1. The 
Wegmann mapping, W , needs to be determined next. 
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Let c (i = 1.2) denote the doubly connected region with 
smooth closed boundaries. Here, the subscript « refer, 
to the inner (i = 1) snd outer [i = 2) boundaries of the 
given jet flow cross-section. It is assumed that t, can e 
parameterized by 

*,.(<) on 0 <«<A ^ » = 1 * 2 < 2 ' 23 ' 


where r. is a smooth regular curve. The aim of the 
method is to determine a mapping g< such that the 
boundaries of the annulus C are mapped onto c.. This 
involves determining the image of a point on the unit an- 
nulus on the boundary of the jet cross-section. As both 
the jet cross-section and the annulus are known smooth 
functions, only the angle 6 of a point on the annulus 
needs to be known to determine the angle, say r(0), of 
its image satisfying: 


( 2 - 24 ) 


The Wegmann method 14 solves for the function, r, 
which is also known as the boundary correspondence 
function. Any real function r(0) satisfying eqn. (2.24) 
and - 

t(6) - £*, fii 6 R (2-25) 

that is 27 r periodic is defined to be a boundary corre- 
spondence function. The image r(0) is determined by 
satisfying eqns. (2.24)- (2.25) . This is obtained itera- 
tively starting from a good guess t and determining the 
correction factor r? such that, 


(2.26) 


For the doubly connected region, there will be two 
boundary correspondence functions for the two edges 
of the shear layer. The interior radius of the annulus 
is also unknown, and must be determined. It is found 
that good initial guesses for both the boundary corre- 
spondence functions and the interior radius are required 
to obtain quadratic convergence. The details of the it- 
erative scheme can be found in Wegmann 14 and Baty 16 . 
The conformal mapping is generated numerically using 
the Cauchy Integral Theorem after the boundary corre- 
spondence problem is solved. 

Once the mapping is completed it is Decessary to 
discretize the governing equations given by (2. 17)- (2.19) 
in the computational domain. A hybrid technique which 
combines a series approximation with a finite-difference 
technique is used here. Let the functions to be approxi- 
mated, in the present case the perturbations associated 
with the shock cell structure, be represented in a series 
of the form: 

4 >{ yi,y 2 ) « XI«(yi.-.y2)/<(yi). ( 2 - 27 ) 

> = 0 


with the coefficient. a(yw,¥2) taken to be function, of 
y, and the basis functions, /, represented m terms of 
Chebyshev polynomials deflned by, 


/.(Vi) 


iii-v?)r w (y,)(-ir i l 

c,N 7 ( yi - yw) 


(2.28) 


r„(y,) is the derivative of the N th order Chebyshev 
polynomial and the constants c, are given by, 

c 0 « c* = 2 and c, = 1 otherwise. (2.29) 

The grid in the yj coordinate direction is deflned by, 


= cos S for 3 = 0,1,2 ,N 


(2.30) 


The model equations can be discretiied readily as the 
basis function, evaluated at the grid points, satisfies re- 
lation . , 

/.( y.y) = 6a ( 2 ‘ 31) 

where % is the Kronecker delta. In the present model 
the derivatives of the approximating basis functions at 
the grid points must also be determined. These deriva- 
tives can be obtained using the relations given by Got- 
tlieb et al. u 

The flow variables are all approximated by a series 
of the form given by (2.27). The approximating senes 
are substituted into the governing eqns. (2.17}-(2.19j. 
These equations are then evaluated at the interior gn 
points. With the use of eqn. (2.31) a system of linear 
ordinary differential equations in the unknown senes c<^ 
efficients is produced. The full form of the equations is 
given by Bhat 16 . 

3. Calculation Procedure 
When the jet shear layer is of finite extent a numeri- 
cal solution to the governing equations must be obtained 
in that region. In order to obtain the numerical solu- 
tion, the solution of the governing equations in regions 
of constant mean flow properties must be found first. 
The separable solutions for the flow variables in a polar 
coordinate system may be found using the technique de- 
scribed by Morris 16 . The solutions for the flow variables 
axe found to be of the form: 


-1 g\ _ _x 2^£l — L—J (jX‘r ) exp(ing) 

P(G*)- jU An lRe + iaM]) K ’ 

(3.1) 

u,M)= £ |il,J,(iVr)+B.J.W|KPM) 

n=_ “ (3.2) 
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Mr.*)- t 

ns-co 

- fl„Y- / " + i (ar) " C«7^(»Ar)]exp(»ntf) 


(3.3) 


u,(r,0) = £ !A n ^ J„(»A*r) + J3 n - Jn+ i(»*0 


+ C„ — (J n (iAr))]exp(«n?) 
dr 


(3.4) 


nere, 


nd, 


(o 3 Re + iaX^hfJ)_ (3.5) 

(Re + t oMy) 


A 3 = o 3 + toRe 


(3.6) 


•uter region: 


p(r,<?)= £ (3.7) 


*.M) = £ ID-rJOiar) g) 

n= — oo 

+ £ n i7i 1) ( i ’ ar )] ex p( tn ^) 

■lr(r,0)= £ [^{^i+\(«r)-r^|jy‘ViWl) 


n = — oo 


(3.9) 


- £'„^'lVi( 1Qr ) “ ^ L ^ 1) ( lQr )]«P(' nS ) 

n= — oo 

+ iE„Hj,+ ] (tor) + — ^-{ifi 1 ) (‘ ar )}) ex P( ,ne ) 

(3.10) 

•ere, (u,, u r> u 0 ) are the components of the velocity Alle- 
gations in the (z, r, 6) coordinates. J n ^ the Bessel 
motion of the Erst kind of order n, is the Han- 

el function of the first kind of order n. These solutions 
^rve as the initial conditions to the governing equations, 
he boundary conditions along the two bounding radial 
nes of the shear layer must also be specified. These 
editions depend on whether the solutions sought are 
id or even about the boundaries. In the present case, 
it a uniform pressure perturbation at the jet exit, so- 
itions are sought that are even about both the minor 
.id major jet axes. The boundary and initial condi- 
ons are applied in the rectangular computational do 
min. The initial conditions are satisfied on the upper 


„d lower boundaries of the recUngle, which correspond 
to the inner and outer edge, of the shear layer The ini- 
tial condition, in (r.f) coordinate, are Uniformed ui 
term, of the computational coordinate, using t e me 
ric generated by the conformal mapping. The boundary 
condition, on the other edge, of the computational dc^ 
main, which represent the two bounding radial lines, are 
.atisfied by .olving for the first and last coefficients of 
the approximating series or its derivative. 

The system of differential equations, derived in the 
last section, may be written as six first- order coup e or- 
dinary differential equation, in the unknown coefficients. 
These equations are integrated from the boundaries for 
each value of n in the series (3 . l)-( 3 .10) . A linear su- 
perposition of these solutions is matched at »ome in- 
termediate location in the computational domain T is 
results in a system of homogeneous equations in the un- 
known series coefficients. The axial wavenumber, are 
determined by teroing the determinant of the coefficien 
matrix using a Newton-Raph.on iterative technique. 

The mean velocity profile is assumed to be given 

W( r ) = Wyexp|-ln(2)rj 3 ] (3- 11 ) 

where r, = n/b, where n is measured normal to the edge 
of the potential core and b is the local half-width of the 
jet mixing layer. For a circular jet, n - r an e 
relationship between the potential core radius h and b 
is obtained from the condition of conservation^ axial 
momentum. This is given by Tam and Moms . In the 
case of the elliptic jet, these values are obtained from the 
measured mean velocity profiles 31 for an elliptic noisle 
of aspect ratio 2 operating at its design condition of 
M rf =1.52. These data include the mean velocity profiles 
along the major- and minor-axes for several downstream 
locations. By fitting the half-Gaussian velocity profile, 
given by eqn. (3.11), the location of the potential core, 
hi (along the major-axis) and h, (along the minor-axis), 
the half-velocity point, tj (along the major axis) and 
fc 3 (along the minor axis) are obtained at the various 
downstream locations. The experimental values of the 
half- widths and the potential core radii are shown in 

Figure 2. . _ 

The half-widths and the potential core radii are 

taken to be a linear function of the axial distance. In the 
calculations for the circular jet, the growth rate of the 
shear layer is taken from the data of Birch and Eggers 
which is 1.266 times the inverse of the spreading rate o 
the mixing layer, c . 

The mean density p may be calculated from the 
mean velocity using a Crocco’s relationship The tur- 
bulent Reynolds number, based on b, the half-width, is 
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where M, and Mt we the fully expanded jet and de- 
,ign Mach numbers. The relationship between the fully 
expanded and design conditions are given by Tam and 
Tanna*. The fully expanded jet is represented by a col- 
umn of uniform flow bounded by a vortex sheet. 

The linearised equation for the pressure perturba- 
tion inside the vortex sheet is, 

vV-M;^=0 (3.14) 

Using elliptic cylindrical coordinates, related to Carte- 
sian coordinates by 

z = a cosh(p) cos(0 ) 

y — CL sin h(p) sin(0) 


Fig. 2 Variation of Shear Layer Parameters with 
Axial Distance of a 2:1 Elliptic Jet. / r y > ^>^2/ r ;j 
o.Aj/ry; My = M d ~ 1.52. 

set equal to 300. The numerical integration is performed 
using a fixed step size fourth-order Runge-Kutta scheme 
with 32 steps in the domain of integration. Calculations 
are performed in one quadrant only based on the as- 
sumed symmetry of the mean velocity profile about the 
major- and minor-axes. In the case of the circular jet, 
three interior lines are used, while for the elliptic jet, 
nine interior collocation points are used. The pertur- 
bation pressure as a function of downstream distance is 
calculated from the pressure perturbation at the nozzle 
exit and the axial variation of the wavenumber, a, for a 
given mode. The pressure may be written as 

(Pru-)o ex P| J ( J ‘ 12 ) 

rvsr0r=3 

where (p n r)o is the amplitude of the n-r mode at the 
nozzle exit and (a nr ) is the complex axial wavenumber 
for that mode. These amplitudes are calculated from a 
vortex- sheet model of the jet. 

The vortex-sheet model was proposed by Prandtl 23 
and Pack 24 , and extended by Tam and Tanna 8 . The 
weak shock celJ structure is modelled as a small- 
amplitude disturbance superimposed on an otherwise 
perfectly expanded jet. The assumption of weak shock 
cells restricts the analysis to supersonic jets operating 
at slightly off-design conditions given by 

| U 7 - M 7 d \ < 1.0 (3.13) 


2 -S 


d 7 p 


eqn. (3,14) becomes 

2 ,aV 

a 2 (cosh 2p — cos 20) dp 7 d& 7 
The vortex sheet is bounded by the ellipse p » Po- 
The boundary conditions are 

D s = 0 on the ellipse p = Po (3-16) 


) + a 2 (M?-l)p S = 0 
(3.15) 


and, 


- — o rj S = Ap everywhere inside the ellipse p — po 
' P y (3.17) 

A separable solution is sought in the form, 

p S (M)»W?(f) (3-18) 

The equations for F(p) and G(f) are found to be, 

— (A — 1q cosh 2p).F = 0 (3-19) 

dp 7 

and, 

4. (> — 2q oos28)G = 0 (3.20) 

where A is a separation constant and q = — 

l)a 2 /4- Equations (3.20) and (3.19) are the Mathieu 
and the Modified Mathieu equations respectively with 
parameters A and q. In general, the solutions of these 
equations are given by the four classes of Mathieu and 
Modified Mathieu functions. The only solution which is 
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symmetric with respect to both s and y axe* is given 

by, 


p S (p,0) = Y2 DnCt it ,(p,q)cci „(M) ( 3 - 21 ) 

n = 0 

where D n is a constant to be determined from the initial 
condition, and ce, n and Ce Jn are the Mathieu and the 
Modified Mathieu function* respectively. 

The boundary condition (3.16) requires that 
Ce 7n (po, 9) = 0 8' ve * root8 9nr. In cal- 

culating po, the dimensions of the fully expanded jet are 
used. The pressure perturbation at the nonle exit is 
then given by, 


p s ( P j) = Ap='£l>2 D n r Cc 7n {p, Qnr )«2n(ff, r ) 

(3.22) 


n = 0 r=l 


where, A p is the pressure difference at the nozzle exit, 
calculated using the one-dimensional isentropic rela- 
tions, and is given by eqn. (2.3) in Tam, Jackson and 
Seiner 10 . By means of the orthogonality property of the 
Mathieu functions, the coefficients D nr are found to be, 


D nr =Ap{ J J Cc 7n {p,q n r)ct 7n (e, q nr ) 

\coih2p — C032&]d6dp^ / 
j j J Ctl n (p,q n r)ccl n (6,q n r) 
\cosh2p — cos26]dffdp^ . 


(3.23) 


The integrals are evaluated using the method suggested 
in McLachlan 25 . In the limit, as the aspect ratio tends to 
unity, eqn. (3.23) reduces to eqn. (2.1b) given in Tam, 
Jackson and Seiner 10 . The initial mode amplitude in 
eqn. (3.12) is related to D„ r by, 

(Pnr)o = D nr Ce 7n (p, 9nr)e c 2n(0| 9nr), ( 3 - 2 *) 


for given values of p and 6. 

The calculation procedure may be summarized as 
;he following steps, (l) The pressure perturbation at 
-he jet exit is obtained as a Fourier-Mathieu series ap~ 
oroximation, eqn. (3.22). (2) The axial wavenumber for 
a given mode at a given axial location is obtained from 
a solution of eqns. (2.17)-(2.19) using conformal map- 
ping and a hybrid pseudo-spectral discretization. (3) 
The pressure associated with the shock cell structure 
.s calculated as the superposition of contributions from 
±ach mode from eqn. (3.12). 


1 . Results 

This section present* the calculations for the shock 
cell structure in a circular jet and in an elliptic jet o 
aspect ratio 2. In the case of the elliptic jet, the pre ic- 
tions are compared with the experimental data obtained 
recently at NASA Langley 21 . These measurements are 
for an elliptic jet issuing from a noitle of aspect ratio 2 
and operated at several off-design conditions. The de- 
sign Mach number of this nonle is 1-52. 

Two different operating conditions have been con- 
sidered - an underexpanded jet (fully expanded jet Mach 
number, My = 1.64, design Mach number, U A - 1.52) 
and an overexpanded jet (My = 1.36, Md — 1- )• ^ 

should be noted that in the earlier work of Moms and 
Bhat is , two different operating conditions (underex- 
panded and overexpanded) were considered for the cir- 
cular jet. The predictions obtained for these cases were 
compared with the experimental measurements of N<>- 
rum and Seiner 26 and also with the results obtained with 
the multiple-scales model of Tam, Jackson and Seiner . 
The comparison showed very favorable agreement. In 
this study, emphasis is placed on obtaining predictions 
of the shock cell structure for comparison with the ex- 
perimental data for the elliptic jet. 

Table 1. Initial amplitudes for the various modes of 


elliptic jet, My=l«64, M<j— 1.52. 


Mode No. 

(pnr)o/pn 

Mode No. 

(pnrWPo 

• 01 

0.31506 

21 

0.04636 

02 

-0.12058 

22 

-0.03442 

03 

0.07365 

23 

0.02305 

04 

-0.05301 

24 

-0.01739 

11 

-0.09432 

31 

-0.02636 

12 

0.05190 

32 

0.02611 

13 

-0.03351 

33 

-0.01779 

14 

0.02475 

34 

0.01365 


Table 2. Initial amplitudes for the various modes of 


elliptic jet, My = 1.36, Md — 1*52. 


Mode No. 

(pnr )o/ Po 

Mode No. 

(pnr)o/pc 

01 

-0.33231 

21 

-0.05083 

02 

0.12597 

22 

0.03589 

03 

-0.07692 

23 

-0.02413 

04 

0.05534 

24 

0.01821 

11 

0.10153 

31 

0.02967 

12 

-0.05429 

32 

-0.02703 

13 

0.03506 

33 

0.01859 

14 

-0.02588 

34 

-0.01428 


In order to obtain a reasonable description of the 
perturbation pressure, the number of modes to be con- 
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sidcred in eqn. (3.12) hw to be determined. This can be 
achieved by calculating the amplitudes of the Founer- 
Mathieu series coefficients representing the perturbation 
pressure at the nozzle exit, see eqni. (3.22) and (3.23). 
Henceforth, any given mode is designated by the indices 
nr. For example, mode 01 corresponds to n = 0 and 
r = 1. The initial amplitude of any given mode is a 
function of both n and r as well as the jet operating 
condition. This dependence can be seen in Tables 1 and 
2. The amplitudes presented here are the perturbation 
pressure (p nr )o on the jet centerline normalized by the 
ambient pressure It can be seen that for a given n, 
the amplitude decreases as r increases. The amplitudes 
also decrease as n increases for any given r. It is also 
clear that many modes would have to be considered to 
provide a perfectly uniform exit perturbation pressure. 
The variation of the perturbation pressure at the noz- 
zle exit along the major- and minor-axes is shown in 
Fig. 3 for the underexpanded elliptic jet, where L ^ and 
Lb are the dimensions of the fully-expanded major-and 
minor axes, respectively. Here, contribution from the 
modes n = 0, 1, 2 and 3 are considered where for each 
n, four roots (i.e, r = 1, 2, 3 and 4) are considered. 
This figure reveals that in spite of considering so many 
modes, there is still some nonuniformity in the pres- 
sure variation across the cross-section of the jet. This is 
characteristic of the difficulty of approximating a step 
function with Bessel functions, in the circular jet case, 
and modified Mathieu functions, in the elliptic jet case. 
However, as the amplitude of the higher modes is reason- 
ably small, it is assumed that a practical approximation 
can be obtained by considering fewer modes. Hence, in 
all subsequent calculations for the elliptic jet, the modes 
considered are 01, 02, 03, 04, 11 and 12. The difference 
between the sum of the contribution of these modes and 
the pressure difference at the exit, obtained using isen- 
tropic relations, is less than 10% for both the underex- 
panded jet and the overexpanded jet. Calculations at 
other operating conditions for both circular and elliptic 
jets have shown that the number of modes required to 
obtain a given degree of accuracy is a strong function of 
the jet operating conditions. 

Figure 4 shows the axial pressure distribution ob- 
tained for the underexpanded jet along the centerline of 
the elliptic jet. Figure 5 is for the case of the overex- 
panded elliptic jet. The sum of the contributions from 
the six chosen modes is compared with experimental 
data of ref. 21. Here, the normalized pressure pertur- 
bation is given as a function of downstream distance, z, 
referenced to the equivalent radius of the fully expanoed 
jet, r y As can be seen, in both the cases, there is favor- 
able overall agreement between the measured and calcu- 
lated pressure distributions. The shock cell spacings and 


C . 3 



T , rfie 


Fig 3 Variation of Pressure at the Nozzle Exit for 

Elliptic Jet, My=1.64, M d =l-52. , along the major- 

axis; , along the minor-axis. 


i 



Fig 4 Centerline Axial Pressure Distribution for El- 
liptic Jet, My=1.64, M d =1.52. , present calcula- 

tions; , measured data of ref. 21. 

the pressure amplitudes are in fairly good agreement. It 
should be noted that the numerical results presented in 
Figs. 4 and 5 have been calculated by shifting the results 
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1.0 r- 



Fig. 5 Centerline Axial Pressure Distribution for El- 
liptic Jet, M, = 1.36, M d = 1.52. , present calcula- 

tions; , measured data of ref. 21- 

by about half an equivalent diameter in the axial direc- 
tion. This is done to account for the unknown location 
of virtual origin of the shear layer. Clearly, the present 
model being relatively simple, cannot effectively model 
the initial development of the jet as the flow transitions 
from the nozzle dimensions to the fully expanded jet. It 
should be noted that the potential core is much shorter 
for the elliptic jet compared to a circular jet at the same 
operating conditions. For example, from Fig. 2, it can 
be seen that /i 3 b 0 for z/ry es 12. Also, the calcula- 
tions in Fig. 4 stop at z/ry = 4.5. Beyond that point 
convergence could not be obtained for the eigenvalues of 
the higher order modes. 

In Figs. 6 and 7, the fundamental mode (mode 01) 
for the two cases are compared with the measurements. 
No axial shift has been applied to the predictions in 
Fig. 7. These hgures show that the shock cell spacing 
can be approximated to a reasonable extent by the fun- 
damental mode alone. However, the amplitudes of the 
shock cells cannot be predicted by this mode alone. In 
order to get a better description, a greater number of 
modes as in Figs. 4 and 5 need to be considered. 

As mentioned earlier, non-circular jets have been 
considered with a view to achieving better aeroacoustic 
characteristics. Thus, comparisons of the axial pressure 
distribution of circular and elliptic jets operating under 
identical conditions have been made. Figures 8 and 9 
show this comparison for both the underexpanded and 


0 - 

- 


0.5 - 



*A, 


Fig 6 Centerline Axial Pressure Distribution for El- 
liptic Jet, My=1.64, M d = 1.52. , mode 01; , 

measured data of ref. 21 


!.° p 

r 



-A; 


Fig 7 Centerline Axial Pressure Distribution for El- 
liptic Jet, My=1.36, M„ = 1.52. , mode 01; , 

measured data of ref. 21 

the overexpanded jets. In the case of circular jet, the 
axial pressure distributions are made up of a linear com- 
bination of the first six modes of the shock cell model. 
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These modes are given by the various teroes of the Bessel 
function, J 0 . and thcir amplitudes at the jet exit is given 
by eqn. (2.1) in Tam, Jackson and Seiner 5 * * * * 10 * * . These fig- 
ures show that the shock cell amplitudes as well as the 
shock cell spacings are less for the elliptic jet case. This 
suggests that the amplitude of broadband shock asso- 
ciated noise might be reduced in Lhe elliptic jet case. 
However, since the decrease is relatively small in abso- 
lute terms it is likely to have a negligible effect on the 
radiated noise calculated in decibels. 



Fig. 8 Variation of Axial Pressure Distribution with 
Aspect Ratio, r/d=0.0, My =1.64, ^<*=1.52. , 2.1 

Elliptic Jet; , Circular Jet. 

5. Summary 

In this paper a linear model has been used to predict 
the shock cell structure in non-circular jets. A confor- 
mal mapping technique with a pseudo-spectral hybrid 
scheme has been used to calculate the ■wavelength and 
decay rate of the shock cell modes. Predictions have 
been obtained for a circular jet and an elliptic jet of as- 

pect ratio 2. The numerical scheme used here has been 
shown to be more successful than the earlier scheme of 

Morris and Bhat 13 in obtaining converged solutions for 

the elliptic jet. The axial variation of the various modes 
contributing to the shock cells for the elliptic jet has 

been observed to behave in a fashion similar to those 

of a circular jet. The amplitudes and the shock cell 
spacings for the elliptic jet have been found to be less 

than those for a circular jet for identical operating con- 

ditions. These changes are likely to result in negligible 



Fig. 9 Variation of Axial Pressure Distribution with 
Aspect Ratio, r/d=0.0, My = 1.36, Md = 1.52. > 2.1 

Elliptic Jet; , Circular Jet. 

direct benefits in noise reduction. However, additional 
benefits of the elliptic geometry could result due to a 
modification of the jet’s turbulent structure or a reduc- 
tion in the supersonic region of the jet. These mecha- 
nisms are being explored by the authors. 
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Abstract 

This paper presents a detailed account of the hydrodynamic stability character- 
istics of the initial region of an elliptic jet. A realistic mean velocity profile is used. 
Calculations of growth rates, phase velocities and eigenfunctions are presented. The 
growth rates of all modes in the initial mixing region are found to depend on the 
minimum momentum thickness. Pressure fluctuations are found to be greatest for 
all modes close to the major axis. An irregular normal mode is found at larger 
eccentricities. All modes, odd or even about the major axis and with periods of tt 
or 27 r have similar growth rates in the initial mixing region. 
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Introduction 


This paper presents a detailed account of the hydrodynamic stability character- 
istics of an elliptic jet in the initial mixing region. The results of the analysis serve 
two purposes. First, they provide a reference case for the verification of analyses 
that consider jets of arbitrary shape. Second, they provide some insight into the 
initial development of large-scale coherent structures in a turbulent, elliptic jet. In 
a previous paper, ref. 1, the author presented some preliminary calculations. The 
present paper extends those calculations and corrects some misinterpretations. 

Purely round jet geometries are the exception rather than the rule in practical 
applications. Jet engine exhausts are fitted with mixing devices to reduce noise 
and decrease the length of the exhaust plume. Non-circular exhausts also occur 
in V/STOL applications in enhanced lift and thrust-vectoring devices. The en- 
hanced mixing properties of non-circular jets make them attractive components for 
fuel-injection and high-speed combustion. A number of recent experimental investi- 
gations have been conducted to examine the properties of non-circular jets 2 8 . At 
low Reynolds number an elliptic jet develops in an unusual fashion as the major 
and minor axes of the jet switch several times with downstream distance. This may 
be associated with the mutual interaction of adjacent elliptic toroidal vortices. At 
high Reynolds number the number and location of any axis-switching remains un- 
clear. However, in this case, experience has shown that the gross properties of the 
large-scale structures in the turbulent mixing layer may be modeled as instability 
waves. This approach has led to a better understanding of jet mixing noise radia- 
tion in high speed jets 9 , and the effect of acoustic excitation on the development 
of turbulent jets 10 . Detailed comparisons between the predictions of instability- 
wave models and turbulence in mixing layers 11 and wakes 12 have also been made. 
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Thus a knowledge of the stability characteristics of an elliptic jet should pro\ ide 
description of the gross features of the large-scale structures in such a flow. 

Crighton 13 examined the stability of an elliptic jet with a “top-hat” velocity 
profile. He obtained some solutions for large eccentricity, in which limit the jet had 
stability characteristics similar to those of a two-dimensional jet. Calculations for 
a wide range of eccentricities for this vortex sheet representation were obtained by 
Morris and Miller 1 . Their numerical results supported Crighton’s asymptotic solu- 
tions. However, the vortex sheet approximation is only valid in the low-wavenumber 
limit. It indicates instability at all frequencies. The observed limited bandwidth of 
unstable frequencies is determined by the finite width of the mixing region. Only 
by using realistic mean velocity profiles can a most-unstable or neutral frequency 
be obtained. Thus in the present paper the stability of an elliptic jet represented 
by a continuous axial velocity profile will be considered. 

A numerical method for the calculation of the stability characteristics of jets 
of arbitrary shape has been developed by Koshigoe and Tubis 14 . Their calculations 
for an elliptic jet compared favorably with the earlier calculations of ref. 1. The 
integral approach that they use suffers from limitations in accuracy as will be shown 
below. However, the technique does not depend on the separability of the stability 
equations and is thus a very attractive approach. 

In this paper the stability characteristics of elliptic jets are documented for 
several eccentricities. The numerical solution requires the evaluation of modified 
Mathieu functions for arbitrary complex argument. Both eigenvalues and eigenfunc- 
tions are presented. In the next section the stability equation in elliptic cylindrical 
coordinates and the asymptotic forms of solution that satisfy the boundary condi- 
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tions are developed. The numerical evaluation of these asymptotic solutions is then 
described. The stability characteristics of an elliptic jet in the initial region are then 
given. Finally, the relationship between these calculations and the development of 
a realistic elliptic jet is discussed. 


Analysis 

A jet flow is considered issuing from an elliptic nozzle. The problem will be 
developed in elliptic cylindrical coordinates (p,0,z). These are related to the 
Cartesian system by: 

x = a cosh p cos 6 

y = a sinh p sin 6 (1) 

z — z . 

The jet axis is aligned with the & direction and the axial mean velocity of the 
jet W[p,6) is assumed to be a function of p and 6 only. This is the parallel flow 
approximation of hydrodynamic stability. This assumption leads to the leading 
order problem in a multiple-scales analysis to include the effects of flow divergence. 
These effects are not considered in the present analysis. 

A Poisson equation for the pressure is obtained by taking the divergence of the 
momentum equation and using the equation of continuity. The resulting equation is 
linearized about the mean flow. The velocity fluctuations are eliminated in favor of 
the pressure fluctuation using the linearized momentum equations. If the pressure 
fluctuation is written in the form, 


p{p,6,z,t) = p{p,6)exp\i{az - wt)], 


( 2 ) 
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then the equation for p may be written, 


a 2 p &_p 2 a(dW_dp dW dp 

dp 2 dO 2 H \ dp dp dd d6 


~~ jcosh(2p) — cos(20)] p — 0, (3) 


where H = u) — cxW . 


If W is taken to be a function of the “radial” coordinate p only then a separable 
form of solution may be sought. This separation is valuable in that highly accurate 
stability calculations may be performed without excessive computation. However, 
as will be seen below, this results in a link between the momentum thickness distri- 
bution of the mean flow' and the eccentricity of the jet which may not be physically 

realistic. 

If a solution for p is sought in the form p = J?(p)T(0) then T and R are found 
to satisfy, 

^ + |A-2,cos(2«)]T = 0, W 

d? R 2 or d,W dR , . P 0 (5) 

and — + — — -j- - - 2 <? cosh(2p)J/t - 0, 

dp 2 n dp dp 

where q = -a 2 a 2 /4 and A is a separation constant. In general q is complex. 
Equation (4) is Mathieu’s equation. Equation (5) is the Rayleigh equation in elliptic 
cylindrical coordinates and reduces to the modified Mathieu equation in regions 
where W is constant. 


The solutions of eqn. (4) are of four types that are odd or even in 6 and w ith 
period n or 27T. Details of the evaluation of the Mathieu and modified Mathieu 
functions and their characteristic numbers are given in refs. 15 and 16. The no 
tation given by Abramovitz and Stegun 16 is used in the present analysis. The 
characteristic numbers were obtained numerically as the eigenvalues of the matrix 
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for the coefficients of the sine or cosine series of the Mathieu functions (Bee ref. 16, 
eqns. 20.2.2-20.2.11). The series was truncated after 11 terms. This gave 8 dec- 
imal places of accuracy for all the values that could be compared with tabulated 
values 16 . The modified Mathieu functions were obtained from a series of products 
of Bessel functions (see ref. 16, eqns. 20.6.7-20.6.10). 


Only certain combinations of these functions permit the physical boundary 
conditions to be satisfied. These were determined by Crighton 13 . In the ambient 
fluid surrounding the jet the pressure fluctuation must vanish as p 1 oo. This leads 


to the following forms of solution. 


p{pj) -» 


A ce 2r+P (*) M4 3 r } +p (p) 

A se 2r + r (0) Ms £!,». 


( 6 ) 


Mc ( 3) and Afs (3) are Mathieu-Hankel functions. If p = 0 the solutions are of period 


7 r , if p = 1 the solutions are of period 27 t. 


The interfocal line p = 0 extends from (x,y) = (-a,0) to (a,0). If the pressure 
and velocity components are to be continuous across the interfocal line then the 
solution for p(p,$) for small p must take the form, 


-I 0\ f B Ce 2r + p(0) Ce 2r + p (p) 

p(p ’* ) = \ Bse ir + P (B) S e 2r+P (p). 


( 7 ) 


This gives the asymptotic form of the solution within the potential core of the jet. 

Crighton 13 considered the case of a vortex-sheet representation of the jet flow. 
That is, 

_ I for » < (8) 

’ [p > \ 0 for P > PC- 

Continuity of pressure and particle displacement at the vortex sheet requires that. 


a Ip! = o 


dp 
d P . 


= 0, 


( 9 ) 
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where A| ) denotes the change across the discontinuity. Application of these condi- 
tions for the forms of solution given by eqns. (6) and (7) leads to a set of dispersion 
relationships. For the even modes these take the form, 


C4+>o)M‘£ ) +>o) „ [, 

Ce 2 r + p(Po)Mc 1 


( 10 ) 


where primes denote differentiation with respect to p. A similar result ma) 
tained for the odd solutions. Crighton 13 obtained asymptotic solutions for these re- 
lationships for large eccentricities. Morris and Miller 1 solved eqn. (10) numerically 
for a wide range of eccentricities. They showed that as the eccentricity increased the 
growth rates of the even modes decreased but those of the odd modes, in particular 
2r + p = 1, increased. This could be interpreted to mean that as the jet's eccentric- 
ity increases the preferred mode would switch from being in phase around the jet, 
such as the axisymmetric mode in the round jet, to being out of phase about the 
major axis, such as the antisymmetric mode in the two-dimensional jet. However, 
these conclusions should only be valid at very low frequencies in a real jet. These 
deficiencies are addressed below where a mean velocity profile with finite thickness 
is considered. 


Calculations 


The mean velocity profile 


The mean velocity profile considered is analogous to that chosen by Michalke 
to describe the initial mixing region of a circular jet. For the elliptic jet the profile 
is taken to be 


W{p) 


1 o < p < p x (11) 

\ |1 + tanh {B(l - sinh p/ sinh p o )/20b}] f° r P > P< 
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A and B are the lengths of the semi-major and semi-minor axes of the ellipse 
defined by p - po and 8 B is the momentum thickness along the minor axis given 

by, 

r °o 

6 B = W(1 - W)dy : x = 0. I 12 ) 

Jo 

In eqn. (12) p, is chosen such that tanh {B{\ - sinh p./sinh p 0 )l28 B ) > s close to 
unity and p 0 is the half-velocity point. Since the minimum valu- of p is zero 
this condition can be met if 6 b /B is sufficiently small. The velocity profile given by 
eqn. (11) reduces to Michalke’s 17 profile in the case of a circular jet. The profile was 
chosen by Michalke on the basis of comparison with experimental data rather than 
being a solution of the equations of motion. As such it is a local representation, 
consistent with the parallel flow approximation of the stability analysis. Axial 
variations are included parametrically through the dependence of 6 B and B on 2 . 
The validity of this choice of profile may be seen by comparing with the experimental 
data of Ho and Gutmark 8 . Figure 1 shows a comparison made at two axial locations 
z/Aq = 0.5 and 2.0, where A 0 is the semi-major axis length at the jet exit. The 
local values of A/B based on the locations of the half velocity points are 1.88 and 
1.49 respectively. The corresponding values of 6 b /B were found to be 0.044 and 
0.223. The agreement between the analytical profile and the experimental data is 
reasonably good except in the inner part of the mixing region along the minor axis 
for the downstream location. 

The momentum thickness along the major axis is defined by 

6 a =[ W{l-W)dx : y = 0. (13) 

Jo 

It can be shown that, 

A6 A = B8 b [l + terms of order {6 b !B) 7 ). (14) 
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Thus as the eccentricity of the jet increases the momentum thickness on the minor 
axis becomes greater than that on the major axis. This is a consequence of requiring 


that IV be a function of p only so that a separable solution can 


be obtained. 


In 


practice this linkage between eccentricity and momentum thickness ratio will not 
exist. In that sense the profile is somewhat unrealistic. However the benefits of 
obtaining highly accurate solutions to the separable problem are felt to justify the 


present approach. 


Numerical method 


The Rayleigh equation (5) was solved numerically using a variable step 
Runge-Kutta algorithm (IMSL routine DVERK modified for complex arithmetic). 
For a given value of frequency and initial guess for the wavenumber two integrations 
were performed starting from „ and n , at which IV = 0.001 and 0.999 respectively, 
toward the center of the mixing region. The starting conditions were based on the 
asymptotic solutions (6) and (7). The resulting numerical solutions at P = Po are 
denoted by R,(p o) and JJ 2 (p o). The two solutions and their derivatives much match 


at p = po, that is, 


Ci-Ri(po) - CiRiiPo) - o 


C\R\ (po) — C'z-^tIpo) — ® 

If eqn. (15) is to have a non-trivial solution for C\ and C 2 , 


(15) 


V (u>, a) = R\ (po) ^ 2 (po) - *;(po)/Mpo) ~ °- 


(16) 


Newton’s method was used to find the zeroes of V(u-,a) 


and hence the eigenvalues 


a. 

Calculations were performed for both odd and even modes about the major axis 
with periods 7 r and 2 tt. The variables were non-dimensionalized with respect to the 


8 



jet exit velocity, the radius of a circular jet of equal exit area sfAB and the uniform 
density. Unless stated otherwise all calculations were for a momentum thickness on 
the major axis of 8 A = 0.02. However, it should be remembered that the momentum 
thickness on the minor axis varies according to eqn. (14). Calculations have been 
performed for three eccentricities: A/B = 1.001, A/B = 2.0 and A/B = 4.0. The 
first case permits comparison of the results with the circular jet case. The areas 
within the contours of p = Po, which are the contours of the half-velocity points, 
were held constant. 

The ceo mode 

The ce 0 mode corresponds to the axisymmetric m = 0 mode in the circular jet 
case. Figure 2 shows the axial growth rates a, for the ce 0 mode as the eccentricity 
changes. The results for A/B = 1.001 are very close to the circular jet results 18 . 
The maximum growth rate decreases as A/B increases. However the frequency for 
the maximum growth remains nearly independent of eccentricity. Table 1 gi'es this 
frequency and the corresponding wavenumber for several modes and values of A/B. 
If the frequencies for maximum growth are non-dimensionalized with respect to the 
major axis momentum thickness 8 A the values for Umax^A are 0.109, 0.113 and 0.113 
f or a/B = 1.001, 2.0 and 4.0 respectively. This suggests that the initial shedding 
frequency is controlled by the minimum momentum thickness, along the major axis 
in the present case, and is nearly independent of eccentricity. This is confirmed by 
the measurements of Husain and Hussain 2 and Gutmark and Ho 3 - 4 . In the former 
experiments the momentum thickness was nearly independent of azimuthal position 

whereas in the latter experiments 8 B - 0.88 A . 

Figure 3 shows the variation of the phase velocity with A/B for the ce 0 mode. 
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For A / B = 1.001 the mode exhibits phase velocities greater than the exit velocity 
that were discussed by Bechert and Pfizenmaier 19 . In general, increasing the eccen- 
tricity makes the phase velocities less dependent on frequency. Table 1 shows that 
the phase velocity for u mai increases with increasing A/B. 

In the round jet, the amplitude of the azimuthal modes are independent of the 
azimuthal location. However the behavior of modes for the elliptic jet is different. 
Firstly, as discussed further below, the modes are not “spinning” but have a fixed 
phase reference to the major or minor axes. Secondly the azimuthal variation is 
determined by the behavior of the Mathieu function that is the solution of eqn. (4). 
As A/B approaches unity q approaches zero, since a goes to zero. Then 
to eqn. (4) are either sines or cosines. However for larger values of A/B, q is a 
complex number with phase determined by the complex wavenumber a. In this 
case the amplitude of the ce 0 mode is no longer independent of azimuthal location. 
Figure 4 shows the azimuthal variation of \T{6)\ for the various eccentricities for 
the maximum amplifying frequency in each case. For A/B = 1.001 the amplitude 
is nearly independent of 6. This is the axisvmmetric mode behavior for the round 
jet. For A/B = 2.0 and 4.0 the amplitude decays rapidly away from the major axis 
and is essentially zero on the minor axis. This would indicate that close to the jet 
exit the pressure and velocity fluctuations associated with the ce 0 mode would be 
greatest near to the major axis. However, because of the form of velocity profile 
chosen this is also the location of the minimum momentum thickness. For uniform 
momentum thickness around the jet exit it is not clear that this behavior would be 
seen. The variation of R{p) with eccentricity is shown in Fig. 5. The distribution 
is plotted along the major axis relative to the half-velocity point and stretched by 
the major axis momentum thickness 6 A . The distributions, plotted in this way. are 
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nearly independent of eccentricity. Thus those motions observed in the round jet 
will be duplicated in the elliptic jet except that the amplitude may vary azimuthallv. 

The ce 2r modes 

In addition to the ce 0 mode there are other even modes of period 7r that are 
unstable in the initial mixing region. Figure 6 shows the axial growth rates for 
these ce 2r modes for r = 0, 1 and 2 for A/B = 1.001. From this figure and 
the numerical values given in Table 1 it is clear that the higher— order modes ha"\e 
similar but smaller growth rates than the ce 0 mode. The azimuthal variation of 
these modes for A/B = 1.001 is given approximately by cos(2r0). The ce 2 mode 
corresponds to the n i. = 2 double helix in the round jet calculations of Mattingly 
and Chang 20 . However it should be emphasized that all the modes in the elliptic 
jet case are phase-locked with respect to the major and minor axes and are not 
spinning modes. This means that any jet that deviates slightly from axisymmetry 
in the mean is likely to exhibit “flapping” rather than spinning modes. 

Koshigoe and Tubis 14 have also calculated the stability of these higher-order 
modes. However the accuracy of their technique, designed for non-separable prob- 
lems, is limited by their two-dimensional grid. In the present calculation the solu- 
tions are analytic in the azimuthal direction and use a variable step-size to achieve 
a prescribed accuracy in the “radial” direction. However, in spite of the relatively 
crude nature of their grid, the results in ref. 14 show all the qualitative character- 
istics of the present results. 

The se] mode 

The modes that reduce to the helical, n = 1 mode in the round jet case are 
the sej and cei modes. They have azimuthal variations that reduce to sin 0 and 
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cos6 respectively. The variation of the axial growth rates for the se, mode for 
various eccentricities is shown in Fig. 7. For A/B = 4.0 the maximum growth rate 
for the se, mode is slightly less than that of the ce 0 mode. The most-amplifymg 
frequencies for the se, mode are also slightly lower than the ce 0 modes. The value 
of ui moI falls from 0.1091 to 0.0704 (when scaled by the momentum thickness 6 A ) as 
A/B changes from 1.001 to 4.0. The phase velocities, shown in Fig. 8 are always less 
than the jet exit velocity. For A/B = 1.001 the results agree closely with the values 
for the n = 1, helical mode in the round jet. At the most amplifying frequency the 
phase velocities are approximately one half the jet velocity (see Table 1). As the 
eccentricity increases the phase velocities become less dependent on frequency. 

For q = 0, the round jet limit, the se, Mathieu function has its maximum 
value at 6 = ±7t/2. However for other values of q this changes. Figure 9 
the azimuthal variation of |T(*)| for various eccentricities. As was found in the ce 0 
mode case the amplitude falls to zero at the minor axis as the eccentricity increases. 
Thus the pressure and velocity fluctuations associated with this mode will also be 
greater close to the major axis (for the present choice of mean velocity profile). 

The se 2r + i modes 

The higher-order odd modes with period 2tt, the se 2r +, modes, exhibit an 
interesting behavior. For real or small q the characteristic numbers of Math.eu’s 
equation are readily classified. However, for complex q pairs of these characteris- 
tic numbers are equal. The location of these branch points has been examined by 
Hunter and Guerrieri 21 . The asymptotic formulae used to calculate the characteris- 
tic numbers changes as one crosses a branch cut that extends radially outward from 
the branch point in the complex plane (see Fig. 4, ref. 21). Across these branch cuts 
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the characteristic numbers exchange their order. This leads to an unusual behav- 
ior for the se 2 r+i modes for higher eccentricities. The characteristic numbers and 
the associated normal modes may be classified initially according to the eigenvalue 
sequence provided by the IMSL routine. Figure 10 shows the complex value of 
c = u>/a for the first four characteristic numbers in the sequence. The eigenvalues 
fall into several groups. There are three continuous sequences associated with the 
second and first, third and second, and fourth and third characteristic numbers. 
There is clearly another sequence with a different behavior in the c-plane that has 
contributions from all the characteristic numbers. This sequence could be thought 
of as indicating an “irregular” mode, though it does not have the same features as 
the irregular mode described by Michalke 22 . In the present case the higher-order 
odd modes have been classified by number if they fall into the smooth sequences 
shown in Fig. 10 or as irregular if they do not. This leads to the axial growth rates 
and phase velocities shown in Figs. 11 and 12 respectively for A/B = 2.0. The max- 
imum growth rate is still associated with the sei mode. All the modes have similar 
phase velocities in the range 0.5 to 0.8 and the dependence on frequency decreases 
as the mode number increases. Though it appears likely that the appearance of the 
irregular mode is associated with the branch points of the characteristic numbers 
of Mathieu’s equation, the mode-switching does not always occur near the branch 
cuts identified by Hunter and Guerrieri 21 . So the appearance of the irregular mode 
remains unexplained. 

The ce 2 r4i modes 

The ce 2r+ i modes are even about the major axis and have a period 27 t. They 
may be identified with a flapping motion about the minor axis in the elliptic jet. 
The cej mode corresponds to the n = 1, helical mode in the round jet. From Table 
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1 it can be seen that the ce, and se, modes have almost identical characteristics, 
approaching the n = 1, helical mode values, as the round jet limit is approached. 
The ce, and se, modes behave very similarly as the eccentricity changes with the 
exception that the frequency of the most-amplifying mode is nearly independent 
of the eccentricity for the ce, mode. In view of the similarity between the two 
modes the growth rates and phase velocities are not shown but the most-amplifying 
frequencies and wavenumbers are given in Table I. 


Discussion 

The present results, subject to the particular mean velocity profile, indicate that 
no particular mode is dominant in the initial mixing region and that the stability 
characteristics are controlled by the minimum initial momentum thickness, which in 
the present case always lies on the major axis. In a given experiment the se 
a particular mode will depend on external influences, such as intended or unintended 
forcing or a feed-back based on the preferred mode of the entire jet flow field. 

Michalke 23 and others have shown how the process of vortex pairing may be 
simulated qualitatively using linear stability theory. Such calculations can describe 
the initial mixing region of a low Reynolds number jet or an artificially excited 
jet at higher Reynolds number. Consider the ce 0 mode which corresponds to the 
asymmetric mode in the round jet case. In the preceding section it was shown 
that the “radial" eigenfunctions for this mode were independent of eccentricity, see 
Fig. 5. Thus the same processes observed in the round jet should occur in the 
elliptic jet. However, in the elliptic jet the amplitude of these motions was seen 
to be dependent on the azimuthal location. For the ce 0 mode the amplitude 


14 



maximum on the major axis and falls rapidly to zero on the minor axis. Thus it 
could be speculated that a vortex roll-up would occur on the major axis but no such 
motions would be seen on the minor axis. Such vortex motions, which resulted in a 
switching of the jet’s major and minor axes, were observed by Husain and Hussain 
and Gutmark and Ho c . However, these processes occur in regions where the jet is 
developing rapidly and the thickness of the mixing layer is increasing. 

The primary purpose of this paper has been to establish the stability character- 
istics of a non-circular jet. The choice of mean velocity profile enabled the stability 
equation to be separated so that highly accurate solutions could be obtained to the 
resulting ordinary differential boundary value problem. The particular choice of 
velocity profile resulted in an initial momentum thickness that varied around the 
jet, being a minimum on the major axis and a maximum on the minor axis. This 
situation could be different in a given experiment as it was in the measurements 
of Husain and Hussain^ and Gutmark and Ho^. However the results given in the 
preceding sections should serve as a test of the accuracy of numerical methods that 
describe the stability of jet flows of arbitrary shape. 
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Mode 

A/B 

wmax 

Q real 

a imag 

ceo 

1.001 

5.4413 

10.1992 

-5.68521 

ceo 

2.0 

5.6578 

10.1355 

-4.5074t 

ceo 

4.0 

5.6518 

9.0558 

-2.5470; 

ce 2 

1.001 

5.4949 

10.3593 

-5.54831 

te 2 

2.0 

4.4171 

8.3818 

-3.1567; 

ce 4 

1.001 

5.6203 

10.7420 

-5.1763; 

sei 

1.001 

5.4537 

10.2381 

-5.6491; 

sej 

2.0 

5.0106 

9.3223 

- 3.6778; 

sei 

4.0 

3.5184 

6.7366 

-2.2132; 

se 3 

2.0 

3.9380 

7.5734 

-2.7856; 

cei 

1.001 

5.4567 

10.2432 

-5.6517; 

cei 

2.0 

5.6577 

10.0272 

-4.5074; 

cei 

4.0 

5.6517 

9.0556 

-2.5470; 


Table I. Frequencies and wavenumbers for maximum rate of growth. 
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Captions 


Fig. 1 Comparison of eqn. ( 1 1 ) with experimental data 8 : z/Aq — 0.5, O , major axis, 

0 , minor axis; z/Aq = 2.0, □ , major axis, A , minor axis. Equation (ll)> 

, A/B = 1.88, 8 b /B = 0.044; , AjB = 1.49, 6 b /B = 0.223. 

Fig. 2 Variation of axial growth rate with frequency for the ce 0 mode. 

, AjB = 1.001; , A/B = 2.0; , AjB = 4.0. 6 A = 0.02. 

Fig. 3 Variation of phase velocity with frequency for the ce 0 mode. For legend see Fig. 2. 

Fig. 4 Azimuthal variation of the amplitudes of the most unstable ceo eigenmodes. For 
legend see Fig. 2. 

Fig. 5 Distributions of most unstable ce 0 eigenmodes along the major axis. For legend 
see Fig. 2. 

Fig. 6 Variation of axial growth rate with frequency for the ce 2r modes. , r = °5 

, r = 1; , r = 2. AjB = 1.001, 6 a = 0.02. 

Fig. 7 Variation of axial growth rate with frequency for the sej mode. For legend see 
Fig. 2. 

Fig. 8 Variation of phase velocity with frequency for the sej mode. For legend see Fig. 2. 

Fig. 9 Azimuthal variation of the amplitudes of the most unstable se) eigenmodes. For 
legend see Fig. 2. 

Fig. 10 Eigenvalues in the c-plane for different characteristic numbers of Mathieu s equa- 
tion. Eigenvalue sequence number. O . 1 : 0.2: A . 3: O.-l. 

Fig. 11 Variation of axial growth rate will) (*- -i'i< nr > (■•* >1" - , i modes. , 

r — 0; , r — 1; , r = 2 ; Q D , irregular mode. AjB - 2.0. 6 A — 0.02. 
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Fig. 12 Variation of phase velocity with frequency for the se 2r+ i modes. For legend see 
Fig. 11. 
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A Linear Shock Cell Model for 
Jets of Arbitrary Exit Geometry 


P. J. Morris, T. R. S. Bhat, and G. Chen 





Abstract 


The shock cell structure of single supersonic non -ideally expanded jets with 
arbitrary exit geometry are studied. Both vortex sheets and realistic mean profiles 
are considered for the jet shear layer. The boundary element method is used to 
predict the shock spacing and screech tones in a vortex sheet model of a single 
jet. This formulation enables the calculations to be performed only on the vortex 
sheet. This permits the efficient and convenient study of complicated jet geome- 
tries. Results are given for circular, elliptic and rectangular jets and the results are 
compared with analysis and experiment. The agreement between the predictions 
and measurements is very good but depends on the assumptions made to predict 
the geometry of the fully-expanded jet. A finite difference technique is used to ex- 
amine the effect of finite mixing layer thickness for a single jet. The finite thickness 
of the mixing layer is found to decrease the shock spacing by approximately twenty 
percent over the length of the jet potential core. 
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1. Introduction 


There have recently been considerable advances in our understanding of su- 
personic jet noise. The results of this work are well summarized by Tam [l]. In 
this paper he shows how the various mechanisms for noise generation in supersonic 
jets may be related to the turbulence in the jet shear layer and its in eraction with 
the nearly-periodic shock structure of the jet. The contribution to the noise ra- 
diation from the turbulence alone, usually referred to as jet mixing noise, may be 
modelled as the radiation of sound by instability waves following the work of Tam 
and Morris [2], and Tam and Burton |3, 4]. The interaction of the instability waves 
with the periodic shock-structure can result in screech tones as shown by Tam |5]. 
The broadband shock noise is associated with the interaction between the random 
turbulence in the jet shear layer and the shock structure. The important features of 
the random turbulence may be modelled as a random superposition of the normal 
modes of the jet shear layer, Tam and Chen [6]. These normal modes are simply 
the solution of the Orr-Sommerfeld equation with mean flow properties given by 
the jet velocity and temperature fields. The success of this model and prediction 
scheme is remarkable in its ability to predict both the near and far pressure fields 

of the jet. 


With the exception of reference 5, Tam’s work has dealt with jets with a simple 
geometry. However recent interest has focussed on the behavior of non-circular jets 
[7-11]. Such jets have enhanced mixing properties. This is important in turbu- 
lent combustion, particularly at high speeds, and for jet exhaust plumes. Nearly 
rectangular nozzles have important uses in STOL applications for lift control and 

vectored thrust. 


2 


In this paper we consider supersonic jets issuing from jet nozzles with arbitrary 
cross-sections. As noted above, Tam and co-workers have shown that both the 
large-scale coherent structures and the shock structure of jets may be modelled 
using a wave analysis. The essential difference between the two cases is that the 
former are travelling waves while the latter may be modelled as stationary waves. 
This means that we are concerned with solving the equations of hydrodynamic 
stability in which the coefficients depending on the mean flow are arbitrary functions 
in a plane normal to the jet axis. Calculations of the stability of non-circular jets 
have also been performed by Koshigoe and Tubis 1 1 2 . 13] and Morris [14]. In the 
former a Green’s function approach was used and in the latter the mean velocity 
profile for the elliptic jet was chosen such that separable solutions could still be 
obtained in elliptic cylindrical coordinates. 

We will assume that the mean flow properties are slowly-varying functions of 
the axial distance. Tam, Jackson and Seiner [15]. showed that the effect of the slow 
axial variation could be accounted for using the method of multiple scales, but that 
the differences in the predicted shock-structure with and without this effect were 
small. Thus we will assume that the mean flow is locally parallel. 

In this paper the shock cell structure and screech tones of a single supersonic 
jet with arbitrary exit geometry are addressed. The jet mixing layer is modelled 
by both a vortex sheet and realistic continuous velocity and density profiles. In 
the former case the area of the vortex sheet cross- section is taken to be that of the 
fully-expanded jet. This is a jet. in which the mass flux is equal to that at the jet exit 
but the pressure has been equalized with the ambient. This is an approximation, 
even for the small pressure differences considered in the present case. However, a 
sample calculation for the two-dimensional jet using linear pressure/turning angle 



relationships showed that the vortex sheet had an average location given by the 
vortex sheet location for the fully- expanded jet. Though the analysis developed 
in the next sections could be applied to segmented exit geometries, such as the 
multitube configuration, the rapid development, of the flow in such cases would 
diminuish the validity of the model. 

It has been found that the boundary element method is a valuable tool for the 
study of jets of arbitrary geometry represented by a vortex sheet. Other numerical 
methods have been employed to consider the cases where the jet mixing layer has 
finite thickness. In the next section the the analytical development is given for both 
representations of the jet. This is followed by the calculations and comparison with 

experiment. 
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2. Analysis 


A typical cross-section of a single jet is shown in Fig. 1. The flow ma) be 
divided into three regions. In regions I and III the mean flow properties are constant. 
These regions correspond to the potential core of the jet and the ambient fluid 
surrounding the jet respectively. Region II represents the annular mixing region in 
which the mean velocity and density of the jet vary. 

As a first approximation a linear shock cell model can be developed in which the 
mixing region of the jet is represented by a vortex sheet. Such a model was used by 
Tam and Tanna [16] in the circular jet case and Tam |5| for arbitrary geometry jets. 
The next approximation accounts for the finite thickness of the mixing layer and 
the variation in the mean flow' properties. To solve these two problems two methods 
have been developed. The vortex-sheet problem is analyzed using the boundary- 
element method. The finite thickness case is solved using a finite-difference solution. 
The analysis for these two techniques is given in this section. 

2.1 Vortex-Sheet Shock Cel) Model. 

The formulation of the problem is the same as that given by Tam [5], and will be 
outlined briefly. However Tam used an eigenfunction expansion method to analyze 
the arbitrary geometry problem. Such an approach is well-suited to problems in 
which the vortex-sheet conforms to a coordinate line in an orthogonal coordinate 
system such as the circular, rectangular or elliptic jet. Tam gave solutions to the first 
two cases and for the elliptic jet in the limiting cases of nearly circular and highly 
eccentric geometries. However such a method is not suitable for more complicated 
geometries. Thus in the present analysis the solution is based on the boundary 
element method. However the eigenfunction expansion method is used in the present 
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paper to find the solutions in the elliptic jet case for any eccentricity. These results 
are used to verify the boundary element calculations. 

Consider a shock cell system in a jet column bounded by a vortex sheet as 
shown in Fig. 2. For convenience, a Cartesian coordinate system centered at the 
nozzle exit with the :r-axis in the direction of the jet centerline will be used. The 
surface of the vortex sheet bounding the fully-expanded jet is given by So(y, 2) — 0* 
There is no disturbance outside the jet. The linearized equations of motion inside 
the vortex sheet arc: 


do 

/»,-V • v + £/y — = 0. 

OX 

(2.1) 

dv 

(2.2) 


V = <LjP- 

(2.3) 


Pj , Uj and ay are the density, velocity and the speed of sound of the fully expanded 
jet. p, p and v are the density, pressure and velocity associated with the linear 
shock cell structure. 

From eqns. (2.1) to (2.3) it is found that the pressure p satisfies the equation 

0, (2.4) 

with p — 0 on the boundary Sc\(y, z) — 0 and, at x = 0, p — A p inside Sn{y.z) = 0 
and Vj_ = 0. vj_ are the velocity fluctuations normal to the jet axis. 

A general solution of the vortex-sheet shock cell boundary value problem can 
be found by writing the pressure fluctuation as: 

7 >(.r,y,z) = <f>{y,z)cos(kx). (2.5) 
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where k is an as yet unknown axial wavenumber. The equation Tor <j>{y,z) may then 
be written 

v\(t> + 0 2 <p = 0. ( 2 - 6 ) 

where 

, ( a 2 a 2 

= {a^ + a? 

and 

a 2 = (A/ 2 - l)k 2 , ( 2 - 7 ) 

with 4> = 0 on So(y,2) = 0. This is an eigenvalue problem with k as the eigenvalue, 
which is to be determined. 

Consider an arbitrary domain as shown in Fig. 3 where the boundary is divided 
into N panels. 6 = £(«.!, a 2 ) are the locations of the node points and a = a(&x,6 2 ) 
are the locations of the mid points of each panel. 




Approximating d<f>/di / by a constant, say g x on each arc F,, we obtain, 


1 A ' f 

= F ("' i 

2 1=1 Jr ■ 


or, 


where, 


N 


,</>(«,) = j{0) ffj. 


y=i 


„(/?) = [ F(ati | 0 do. 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 


In eqns. (2.10) to (2.13) i — 1,2...., N . 

The approximate eigenvalues 0 are obtained from, 

dct[/z,j]/vx/v — (2-11) 

It now remains to set up the elements of the matrix When i = J the integrand 
contains a singularity and so the evaluation of the integral in must be performed 
carefully. We obtain: 


Mtt (/?) — ” \ L>iH { 


T,„v)(2k\ + — 


4 | \ 2 

v H\ X)(0Lx 


Hn 


0L, 




(2.15) 


where Li is the length of the panel on the arc F t , and Hr. and Hi are Struve 
functions of order 0 and 1 respectively, that can be expressed as a series of Bessel 
functions. When t # j the integrand is smooth and Ihe evaluation of the integral 
in fj.ij may be obtained using Simpson’s rule. 

/4 ]) (/?|a t - Cyl) 


i ( Lj 

Fij{0) - a ) 




4 4 


H^ ) {0\a l -a ] \) + H^{0\^-ij + i\) 


(2.16) 
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The eigenvalues are obtained by a local iterative scheme which is given by, 


Ih+i = 0k - 1 ffilh). 


(2.17) 


where, 

f(l3 k ) = Tr|M _1 (^)/ / (^)l- ( 2 - 18 ^ 

Tr denotes the trace of a matrix and a prime denotes the derivative of the matrix 
with respect to 0 k - This is Newton's method written in a convenient form for matrix 
operations. Thus the eigenvalues can be determined starting from an initial guess. 


2.2 Continuous Mean Profiles 


When region II of Fig. 1 is of finite extent a numerical solution must be obtained 
in the mixing region. The general form of separable solution for the fluctuating 
pressure when the mean velocity and density are taken to be independent, locally, 
of the axial distance is given by, 


p(r,6, x,t) = p(r,0) exp \i{kx - w/.)]. 


(2.19) 


where u> is the radian frequency and a polar coordinate system has been introduced. 
The shock-structure may be associated with the zero frequency solutions, i.e. w = 0. 
The pressure fluctuation is found to satisfy the non-separable form of the Rayleigh 


equation [12], 


d 2 p 1 dp 1 d 2 P 

dr 2 ^ r dr r 2 dO 2 


2k 


f dU dp 


(w - kU) 


dr dr 


1 dV dp 

^99 98 


- [ k 2 - M;(c - 


kU) 2 )p = o. 


(2.20) 


For simplicity the flow has been assumed to be isothermal but the variations in 
mean density may be readily included using a Crocco s relationship. 
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For an arbitrarily shaped jet, V is a function of both r and 6. Tims a separable 
solution in r and 6 is not generally available and one resorts to a numerical solution 
in region II, Along the interior edge of this region and in region I, the potential core 
of the jet, the mean velocity is a constant, U } . In the ambient fluid surrounding the 
jet the properties are also uniform. In this case we take the external velocity to be 
zero. Thus on the boundaries between regions 1 and II, and 11 and III, a separable 
form of solution may be obtained and can be written. 

oo 

Pj = ^n^f.(iAyr)exp(tn<?), (2.21) 

ft — — OO 

and, 

/in = jr iJn/f^HiAnrjexplmfl), (2-22) 

n — — oo 

where 

A* = k 2 - My(w - k Uj ) 2 , (2.23) 

and, 

Xl = k 2 -My. (2.24) 

The subscripts j and 0 refer to the solutions in regions I and III respectively. J n is 
the Bessel function of the first kind and order n, and H^ ] is the Ilankel function 
of the first kind and order n. 

For region II the solution may be obtained numerically using a finite differ- 
ence scheme. The region is divided by N radial lines. If the mean flow possesses 
symmetry about any coordinate directions the numerical integration need only be 
performed in a limited sector. For example, if the mean flow is symmetric about 
both the y and 2 axes the numerical integration is performed in the first quadrant. 
This would be the case for the elliptic or rectangular jets. The derivatives with 
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respect to 0 in eqn. (2.20) are approximated by a three-point central difference for- 
mula. Equation (2.20) then gives a system of coupled ordinary differential equations 
for the values of p along the radial lines. The values of p along the bounding radial 
lines depend on whether the modes sought are odd or even about the boundaries. 
In the the former case p is zero on the boundary, in the latter case dp/d6 is set to 
zero. 


The numerical solution is started at t He interior boundary with one term from 
a finite series of the form (2.21). The values of n are chosen to be compatible with 
the modes sought. For example, n = 2m for modes that are periodic in n such as 
the axisymmetric mode in the circular jet case. The solution is repeated for each 
term in the series as the starting condition. A corresponding set of solutions is 
started along each of the radial lines at the exterior boundary. The sets of solutions 
are matched at some intermediate value of radius. For example along any radial 
line at the matching point, 

(2 ' 25) 

n n 

and 

Z A *&I = 'L B '&"' ( 2 - 26 ) 

n n 

where p n / and p n // denote the numerical solutions started at the interior and exte- 
rior boundaries respectively with starting conditions for the n-th term in the series, 
and primes denote the derivative with respect to r. Application of these matching 
conditions along each of the N radial lines gives 2A homogeneous equations for the 
2N unknowns A„ and B n . These equations may be written in matrix form. The 
eigenvalue is obtained by minimizing the determinant of the resulting matrix. 
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3. Calculations 


Calculations have been performed for both the vortex sheet case and the jet 
mixing layer of finite thickness. In the former case the circular, rectangular, and 
elliptic jets have been considered. Comparisons have been made between the pre- 
dictions, analytical results and experimental data where possible. 

3.1 Vortex-Sheet Shock Cell Model 

The shock spacing is given by the wavelength A., associated with the lowest 
eigenvalue k. For a circular jet the analytical result for 0, given by eqn. (2.7), was 
given by Tam and Tanna |1G|. It corresponds to the first zero of the J 0 Bessel 
function. To determine the accuracy of the boundary element method, calculations 
were performed for the circular jet case with various numbers of panels. The results 
are given in Table I. The numerical result is within one percent of the analytic result 
for 20 panels. Thus in all subsequent calculations the number of panels is fixed at 
20 unless noted otherwise. 

Tam and Tanna |16] noted that the dimensions of the vortex sheet to be used 
when comparing with experimental data should be those of the fully-expanded jet 
and not those of the jet nozzle itself. The ratio of the areas of the fully-expanded 
jet and the nozzle may be obtained from one-dimensional isentropic flow formulas. 

That is, , | 



Md 



Ad 

I! 

| 

l + : 

*r-M] 


where Md and M 2 are the nozzle design Mach number and fully expanded jet 
Mach number respectively and A d and Aj are the areas of the nozzle exit and the 
fully-expanded jet respectively. For a circular jet the ratio of the radius of the 
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fully-expanded jet to the nozzle radius is given by the square root of eqn. (3.1). 
For non-circular jets further assumptions must be made. Tam [5] assumed that 
the jet expands or contracts by an equal amount about its perimeter. Approximate 
formulas may then be found for the dimensions of the fully-expanded jet. For the 
cases of the rectangular and elliptic jets these formulas are eqns. (2.42) and (2.43) of 
reference 5. Alternatively the shape of the jet cross-section could be assumed to be 
unchanged in the expansion so that the ratio of the characteristic dimensions of the 
fully-expanded jet to the nozzle is once again given by the square root of eqn. (3.1). 
However, though without any conclusive experimental evidence, the mixing in the 
initial mixing region near the nozzle might be more likely to move the jet to a more 
symmetric form, i.e. a lower eccentricity in the elliptic case. This trend is given by 
Tam’s formulation , at least in the under-expanded case. This formulation has been 
used mostly in the subsequent calculations. However some predictions based on the 
assumption of no change in the jet geometry have also been made. For the elliptic 
jet the formula for the fully-expanded jet dimensions should be written correctly 


as: 


hi 

L,i 


~Ad ~ *J 2 


(3.2) 


where L-j refers to the fully-expanded jet scales and Lj refers to the nozzle exit scale. 
E(s) is the complete elliptic integral of the second kind and e 2 = (l-L 2 /L 2 ), where 


L a and Lb refer to the semi-major and semi-minor axes of the nozzle. 


For the elliptic jet the shock spacing is given by [5], 

Xf — - l)/\/9m • (3-3) 

where </ 0 i is the smallest root of the modified Mathieu function Ce nm (p.o, <7mn), and 
a - y f L 2 a - L 2 b (x o = tanir '(L a /Lb), ( 3 - 4 ) 
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where La and Lb are the lengths of the semi-major and semi-minor axes of the 
fully— expanded iet. Tam [5] gave approximate values for the roots in the limits 
of very small and very large eccentricity, In the present calculations the roots 
have been determined numerically using the Mathieu function routines described 
by Morris [14]. These calculations confirm the asymptotic formulas given by Tam. 
In addition calulations have been performed for an elliptic jet with nozzle axis 
ratio La! L b - 2. This corresponds to the nozzle used in recent experiments at 
NASA Langley [17]. Figure 4 shows the calculated shock spacing as a function 
of fully— expanded jet. Mach number for the 2:1 elliptic, nozzle. The design Mach 
number of the jet is 1.5. The shock spacing is referenced to the semi-minor axis 
of the nozzle. The numerical results (for A^ = 20) agree well with the analytic result 
given by eqn. (3.3) with the zeroes of the Mathieu function evaluated numerically. 
It should be noted that an axis ratio of 2:1 only occurs at the ideally-expanded 
condition. For M j — l.l the fully— expanded axis ratio is 2.22/8 and for Mj 1.8 
the ratio is 1.7753. 

A relationship between the shock spacing in the jet and the screech tone fre- 
quencies has been developed by Tam, Seiner and Mi [18] based on a weakest link 
hypothesis. The screech tone frequency is given by the simple formula, 

r - Ug/Cl , (3.5) 

r 27T (1 + U c /fl-r<,) 1 

where k\ is the smallest wavenumber of the shock cell system, u c is the convection 
velocity of the excited large scale instability waves of the flow and a .00 is the ambient 
speed of sound. The same formula was developed by Powell |19j and Harper-Bourne 
and Fisher [20] using a different, model. In the calculations it has been assumed that 
the convection velocity is 0.7. This assumption was found to be satisfactory in [18] 
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though it may not be generally valid. Though other values have been proposed for 
the convection velocity the present choice is an adequate approximation. It should 
be noted however that the convection or phase velocity for the instability wave that 
achieves a maximum amplitude in the jet shear layer may be calculated from the 
stability analysis (for realistic mean velocity and temperature profiles). Figure 5 
compares the predicted screech tone frequencies from eqn. (3.5) with the measured 
values. The smallest wavenumber was obtained from the boundary element calcu- 
lations shown in Fig. 4. The agreement between the predictions and measurements 
is very good. It is interesting that the experimental data and the predictions ap- 
pear to agree best close to t lie design Mach number. This would suggest that the 
discrepancies between the predictions and the experiments are related to the choice 
of dimensions for the fully-expanded jet. Preliminary measurements [17] of the 
velocity profiles in the Mach 1.5, 2:1 elliptic nozzle at NASA Langley indicate that 
the ratio of the axes, based on the half-velocity points, remains constant for the 
first two potential core lengths. Also shown in Fig. 5 is the predicted screech fre- 
quency where the aspect ratio of the fully-expanded jet is taken to be the same as 
the nozzle. The area of t lie fully-expanded jet is given by the one— dimensional gas 
flow equations. The agreement between the predictions and experiment is better in 
this case particularly at the over-expanded conditions. It is not clear whether this 
assumption would provide better predictions for all exit geometries. 

The circular and elliptic jets have smooth boundaries. As an example of a 
case with sharp corners we consider the rectangular nozzle. The analytic solutions 
in this case are obtained very easily using the eigenfunction expansion approach 
and the formulas for the shock spacing and screech frequencies are given by Tam 
[5] [eqns. (2.44) and (3.2)]. Figure 6 shows a comparison between the numerical 


15 



results for the shock spacing obtained using the boundary element method, Tam s 
analytic results (based on the assumption of large aspect, ratio), and Powell’s [21] 
experimental results. The agreement between the sets of data is quite reasonable. 
For these calculations, in which the aspect ratio of the jets is large, 34 panels were 
used in the boundary element calculations. 


Figure 7 shows a comparison of the predicted screech frequencies using the 
analytic and numerical results and the experiments of Krothapalli et al [22]. As 
expected the agreement between the predictions and experiments is good. Once 
again 34 panels were used in the calculations. 

It should be noted that at the higher pressure ratios strong shocks exist in the 
plume and the geometry of the plume changes with downstream distance. This 
results in the eventual switching of the major and minor axes. The good agreement 
between the predictions of the vortex sheet model and the measurements would 
suggest that the fundamental wavelength of the jet is relatively insensitive to the 
effects of finite mixing layer thickness, jet growth, and non-axial velocity compo- 
nents. This is also indicated by the calculations given in the next section. 

3.2 Continuous Mean Profiles 

As a demonstration of the effect of finite thickness on the shock spacing cal- 
culations were performed for a circular jet. Clearly this case could be examined by 
using the separable form of solution and integrating along a single radial line. In 
fact this calculation has been performed to test the numerical accuracy of the gen- 
eral scheme. However the purpose of the example was to test the general approach 
and calculations for other geometries must await mean flow data at high speeds. 


16 


The calculations were performed with three evenly spaced radial lines in the 
first quadrant. The mean velocity profile was taken to be, 

V(r) --- Uj exp[- ln(2) ?;*], (3*6) 

where t) = (r — h)/h, h is the potential core radius and 6 is the local half-width of 
the jet mixing layer. The relationship between h and b may be obtained from an 
integral form of the axial momentum equation (see for example, lam and Morris 
[23]). The jet was assumed to be isothermal. The integration of the system of 
ordinary differential equations in the radial direction was performed using a fixed 
step size fourth-order Runge-Kutta scheme with only 20 steps in the integration 
region, 0 < r? < 3.0. This relatively small number of steps appeared to be adequate 
when compared with the results using a variable step-size scheme that used several 
hundred steps. 

Figure 8 shows the variation of the shock spacing, based on the lowest wavenum- 
ber, as a function of the mixing layer thickness. The jet Mach number My = 1.4. 
The wavelength is relatively independent of the mixing layer thickness though the 
shock spacing does decrease as the thickness increases. For small values of b the 
spacing approaches that given by the vortex sheet model. Only the lowset wavenum- 
ber associated with the axisymmetric mode has been sought. Higher axisymmetric 
modes exist corresponding to the zeroes of the Jn Bessel function in the limit of 
the vortex— sheet. Tam, Jackson and Seiner [15] showed that when the contribu- 
tions from all these modes are superimposed a remarkably detailed prediction of 
the pressure distribution, associated with the shock structure, may be obtained. A 
similar result could be obtained with the present method for other jet geometries. 
However, it should be noted that Tam et al’s calculations included the dissipative 
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effect of the shear layer turbulence through an eddy viscosity. This results in an 
axial decay of the strength of the shocks. No such decay would be obtained with 
the present formalism. However the inclusion of the viscous effects would simply 
convert the equation satisfied by the fluctuations from the Rayleigh equation to an 
Orr-Sommerfeld equation and the numerical analysis would proceed in a similar 
fashion. 

The method developed in this section is being applied by the authors to other 
more complicated geometries than the circular jet. There are several modifications 
to the numerical procedure which make the calculations more efficient in jets with 
high aspect ratios. However the basic principle of matching the linearly independent 
solutions along radial lines is the same as described here. These calculations will be 
presented later. The present calculations have served to validate the basic technique 
and also to demonstrate the magnitude of the effect, of finite mixing region thickness. 
It can be seen from the present calculations that the vortex-sheet approach provides 
a remarkably good prediction of the shock spacing when compared with the more 
realistic results for the finite thickness mixing layer. 
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4 . Summary 


This paper has given t he formulations for the calculation of the shock cell spac- 
ing of single jets of arbitrary geometry. Methods have been given for both vortex 
sheet models of the jet shear layer and for continuous mean profiles. The vortex 
sheet calculations for the single jet were performed using the boundary element 
method. This enabled the shock structure to be predicted in jets of arbitrary ge- 
ometry. Calculations have been given for circular, rectangular and elliptic jet cases. 
The effect of finite mixing layer thickness on the shock spacing has been predicted 
using a finite-difference solution in the mixing layer. The shock wavelength has 
been found to reduce as the mixing layer thickness increases. 
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Analytical value of )3 ~ 2.40482 


Number of panels 

/bminenen 

8 

2. 5860 

12 

2.4765 

20 

2.1287 

30 

2.4151 

40 

2.4100 


Table I. Boundary Element Calculations for a Circular Jet. 


Figure Captions 


Fig. 1. Sketch of Regions of .let Cross Section. 


Fig. 2. Sketch of Vortex- Sheet Shock Cell Model. 


Fig. 3. Sketch of Boundary Divided into Elements. 


Fig. 4. Variation of Shock Spacing with Fully-Expanded Jet Mach Number for 
Elliptic Jet. L a jLb — 2.0. , present calculations, N — 20;* , analytic soution, 

[5]- 


Fig. 5. Variation of Screech Tone Frequency with Fully-Expanded Jet Mach Number 

for Elliptic Jet. L a / L\, ~ 2.0. , present calculations, N = 20; * , analytic 

solution, [5]; A , experimental values, [17]; , present calculations assuming 

constant jet aspect ratio. 


Fig. 6. Variation of Shock Spacing with Fully-Expanded Jet Mach Number for 

Rectangular Jet. , present calculations, b/h = 5.83, N = 34; * , analytic 

solution, [5]; A , experimental values, b/h = 5.83, (19). 


Fig. 7. Variation of Screech lone Frequency with Fully-Expanded Jet Mach Number 

for Rectangular Jet . , b/h = 50/5, , b/h — 50/3, present calculations, 

N = 34; * , analytic solution, |5); experimental results, |20], A, b/h = 50/3, Q, 
b/h = 50/5. 


Fig. 8. Variation of Shock Spacing with Mixing Layer Thickness for Circular Jet. 
Mj = 1.4. 




Fig. 1. Sketch of Regions of Jet Cross Section. 
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Fig. 6. Variation of Shock Spacing with Fully-Expanded Jet Mach Number for 

Rectangular Jet. , present calculations, b/h = 5.83; * , analytic solution, [5|; 

A , experimental values, b/h = 5.83, (19). 
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INSTABILITY WAVES IN TWIN SUPERSONIC JETS 


Philip J. Morris 
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University Park, PA 16802 USA 

Calculations are presented for the characteristics of instability waves in the 
initial mixing region of twin circular supersonic jets. Two models for the basic 
jet flow are used. In the first, the jets are modeled as two circular vortex sheets. 
In the second, realistic velocity and density profiles are used. It is shown that 
the unsteady flow fields of the two jets do interact before the time-averaged jets 
flows have merged. The normal modes or instability waves are classified by their 
symmetry properties in the twin jet case and their asymptotic behavior for large jet 
separations. Calculations of the growth rates and phase velocities are made for these 
modes as a function of jet separation and mixing layer thickness. The associated 
pressure distributions are also presented. In the realistic jet profile calculations the 
effect of jet separation is found to be relatively weak. For modes that are even 
about the symmetry plane between the two jets the pressure levels are found to 
increase near this plane as the jet separation decreases. 
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1.0 Introduction 


When supersonic jets from convergent-divergent nozzles operate at off-design 
conditions they can produce intense screech tones. Powell (1953) made early ob- 
servations of this phenomenon and proposed a feedback mechanism for the screech 
tone production. More recent experiments and analysis by Tam, Seiner and Yu 
(1986) showed that the feedback loop consists of downstream propagating large 
scale structures in the jet mixing layer that interact with the shock cell structure 
to generate upstream travelling acoustic waves. If these acoustic waves trigger ad- 
ditional flow disturbances at the jet lip with the correct phase then the feedback 
loop is established. Analyses, based on this model, have made excellent predictions 
of screech tone frequencies in both circular and non-circular jets; see Tam (1986) 
and Morris, Bhat and Chen (1989). 

Seiner, Manning and Ponton (1988) showed experimentally that for two closely- 
spaced supersonic jets, operating off— design, the dynamic loads associated with 
the screech tone can reach levels, upstream of the jets’ exits that could result in 
structural damage. Tam and Seiner (1987) noted that the screech tone frequency 
of the twin jets was slightly different to that of the single jet and that the acoustic 
intensity in the inter-nozzle region exceeded that of the direct sum of two non- 
interacting screeching jets. This suggests that there is a strong interaction between 
the unsteady flow and acoustic fields of the two jets. The analysis and calculations 
described in this paper help to quantify the effects of jet separation and operating 
conditions on the nature of this interaction. 

Turbulent mixing in free shear flows is controlled by the dynamics of large scale 
coherent structures. The local characteristics of these structures may be described 
by linear instability theory. This has been demonstrated by the experiments of 
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Gaster, Kit and Wygnanski (1985), and Petersen and Samet (1988) among others. 
In their experiments they compared predictions of the amplitude and phase of the 
axial velocity fluctuations, based on linear stability theory, with phase-averaged 
measurements in an excited shear layer and a jet. The agreement between predic- 
tions and experiment was very good though only the local distributions and not the 
amplitude were predicted. This close agreement between the predictions of linear 
stability theory and the properties of the large-scale coherent structures has formed 
the basis for theories of turbulent mixing and supersonic jet noise. For example, 
Tam and Morris (1979), and Tam and Burton (1984a, b) predicted the noise radi- 
ation from instability waves in supersonic shear layers and jets and obtained very 
good agreement with experiment. 

For supersonic jets the three main components of noise radiation are, turbulent 
mixing noise, broadband shock associated noise and screech. In each case, the 
essential component of the turbulence responsible for noise generation are the large 
scale structures. It should be noted that this is not the case for subsonic jets where 
a complete theory for noise generation and radiation is not available. Tam (1987) 
showed how predictions could be made for each noise component in a circular jet 
using an instability wave model for the large scale structures. 

In the present paper the properties of the instability waves or large scale tur- 
bulent structures in the initial mixing region of twin circular supersonic jets are 
determined. Two models for the basic jet flows are used. In the first, the jets are 
modeled as two circular vortex sheets. In the second, realistic mean velocity and 
density profiles are used. Though the former model fails to provide quantitative 
results it does help to explain the observed modes of instability and interactions 
predicted by the more realistic model. The calculations examine whether the in- 
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stability growth rates, and hence the amplitudes of the large scale structures, are 
modified a s the jet separation and operating conditions vary. In addition, the cor- 
responding changes in the instability wave phase velocity are predicted. It is shown 
that the unsteady flow fields associated with the instability waves do interact be- 
fore the time-averaged jet flows have merged. However, this interaction is relatively 
weak for the operating conditions considered. In Section 2 the general equations of 
motion and analytic solutions common to both models are developed. The details of 
the vortex sheet model and its predictions are then described in Section 3. Section 
4 contains the numerical procedures and calculations for the realistic jet profiles. 
Finally, the role of these predictions in and their relationship to experimental ob- 
servations of the twin plume resonance phenomenon are discussed. 

2. ANALYSIS 

Consider the two circular jets shown in Fig. 1. The time— averaged jet flows are 
assumed to be symmetric about the x — z and x — y planes, where the x-coordinate 
is normal to the jet exit planes. The centers of the jets are separated by a distance 
2h. Throughout this analysis the variables are non-dimensionalized with respect 
to the jet velocity uy, jet density p } and jet radius ay. These values are taken to 
be the fully-expanded jet properties as defined by Tam and Tanna 11 . These values 
are described below. In the annular mixing regions of the two jets, before they have 
merged, there are three flow regions. In region I, the potential cores of the jets, the 
mean velocity and density are constant. Region III represents the stationary fluid 
surrounding the jets. In region II, the annular mixing region, the mean velocity u 
and density p are variable. Polar coordinate systems are introduced and 

(r 2 ,# 2 ) with origins on the jet centerlines. The mean velocity and density of each 
jet are assumed to be a function of their radial coordinates only. This is the locally- 
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parallel flow approximation. The potential cores have radii i?i and the outer edges 
of the mixing regions have radii R ? . The mean static pressure is assumed to be 
constant. The large scale coherent structures are modeled as instability waves. 
Their behavior is governed by the unsteady, linearized, compressible equations of 
motion. Thus, for example, in either polar coordinate system separable solutions 
for the pressure fluctuation are sought in the form: 

p(x, t) = p(r) exp[i(A:x + nO - wt)), (2.1) 


where k is an axial wavenumber, n is an azimuthal mode number, and w is a radian 
frequency. The radial variation of the pressure fluctuation is then found to satisfy 
the equation: 


^ + /I 

dr 2 l r 


1 dp 2k du'i dp r^ 2 

p dr ^ Q dr / dr l 


- pfl 2 Mf + 



(2.2) 


where, 


H = u> — ku, 


and Af 2 = u 2 /c 2 , where Cy is the fully— expanded jet speed of sound. Equation (2.2) 
reduces to Bessel’s equation in regions of constant mean velocity and density. 


A solution for the pressure fluctuations outside the jet mixing layers in region 
III may be obtained in either polar coordinate system in the form 


where, 


p(r,0,x,t) 


^ B n H { n 1] (i\ 0 r) exp[t(fcx -ut + n6)}, 

= — OO 

A 0 = {/c 2 -p 0 u, 2 m;}>. 


(2.3) 


p 0 is the non-dimensional mean density in the ambient medium which is equal to 
the jet static temperature ratio, Tj/T 0 . The branch cuts for A 0 are chosen such 
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that, 


~ 2 * - ar S A ° < 2 *’ 

This ensures that the solutions decay as r — ► oo or are outgoing waves, for positive 
frequency. The symmetry properties of the mean velocity and density field, for 
example, 

u(y,z) = u(-y,z) - u(y,-z) 


indicate that the eigensolutions should be odd or even about the x — y and x — z 
planes. From the latter symmetry property (the former symmetry property is used 
below) the general solution for the pressure fluctuations in the outer region may be 
written, 


OO 

p(ri,r 2 ,6 1 ,6 2> x,t)= £ B n {if' 1 ’ (tA 0 r 2 )e*" 0 ’ 

| exp[i(/tx — wt)], 


n = — oo 
in ( ir - ) 


(2.4) 


±Hl 1 ) (iX 0 r 1 )e i 

where the choice of sign depends on whether the solution is to be odd or even about 
the x — z plane of symmetry. 


It should be noted that this solution does not just represent the sum of the 
contributions from two, non-interacting, individual jets, though the form could then 
be the same. It is simply a convenient form of the separable solution in the outer 
region. The influence of the second jet is included when this general form of outer 
solution is matched with the solution in the interior of each jet. In the present case 
the fluctuations in each jet are affected not only by the outgoing solutions, that 
would exist for an individual, isolated jet, but by the incoming solutions from the 
second jet. These two contributions are included in the outer solution (2.4) 


In order to match the outer solution (2.4) with the pressure fields in regions I 
and II it is convenient to write the solutions in region III in terms of only one of 
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the two polar coordinate systems. This may be accomplished using Graf’s Addition 
Theorem [see, Tranter (1968)]. For example, we may write, 

oo 

ir;*'(a 0 r 1 y'-’ ,, = £ < 2 - 5 ) 

$ = — CO 

Then (2.4) may be written, 

p(r 2 ,0 2 )= ^ B,P, n (iA 0 r 2 ), (2.6) 

n — — oo #= — co 

where, 

/?.„(£) = ± H l n l l{2iXoh)J n {t), (2.7) 

and 6 tn is the Kronecker delta function. It should be noted that 

= (-x)'A.. (2- 8 “) 

and, 

=(-l)*/?.,-„. (2-86) 

The influence of the second jet is seen in (2.7). The first term on the right 

hand side represents the outgoing waves from the jet at y = —h. The second term 
represents the incoming waves from the jet at y — h. For large h this latter term is 
very small and the interaction between the jets is very weak 

The normal modes for the pressure fluctuation given by (2.6) may be separated 
further into modes that are odd or even about the x — y plane. This is accomplished 
by first setting s = -s and n = — n and adding the resulting equation to (2.6) . Then 
the pressure fluctuation may be written, 

OO OO 

p( r 2,07)= [ F * cos(n# 2 ) + iG, sin(n0 2 )J /?.„ (i A 0 ^ 2 ) > (2.9) 
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where, 


F. =\B. + (—!)* B-.J/2, 


(2.10a) 


and, 

G, = {B. -(-l)'B_.]/2, (2.106) 

A similar solution may be found for the pressure fluctuations in the potential core 
region. This may be written, 

OC 

p(r 2 ,0 2 )= [ An cos(n0 2 ) + iC n sin(n0 2 )J J n (»^i r 2 ), (2.11) 

n — — oo 

where, 

X] — k 2 - (u> - k) 2 Mj. (2.12) 

With the form of the solutions known inside and outside the jet the eigenvalues k 
may be determined by matching these solutions at either the vortex sheet location 
or through the finite mixing layer. The former matching is described in the next 
section. 

3. VORTEX SHEET MODEL 


3.1 Analysis 

In this representation of the jet flows the finite mixing layers are replaced by 
cylindrical vortex sheets of unit non-dimensional radius. Across the vortex sheet we 
require continuity of pressure and particle displacement. The matching conditions 
require that, 

= 0 and A [p] = 0, (3.1) 

where A[ ] denotes the change in the argument across the vortex sheet. If the 
interior solutions given by (2.11) are matched with the exterior solutions, (2.9) for 


dp / dr 2 
pQ 2 ~ 
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all n we obtain, for the solutions that are even about the x — y plane, 

y~] F, {py^J'Ai\i)P.n{^o) — {u - fc) 2 A 0 J„(tA 1 )/?' n (iA 0 )| = 0, 

• = oo 

n = — oo, • • • , oo, 

where, 

«.(0 = *..«!“'(£) ± frlV.(2>A.k)-':(e)- ( 3 - 3 ) 

If we use the symmetry properties of given by eqns. (2.8), and note that, 


F. =(— 1 )*F_„ 


(3.4) 


then the independent equations yielded by (3.2) may be written, 

f>. {r n <5.„ ± [jsriV. (2*Ao/t) + (-l)*^ 1 _ > .( 2 tA 0 / l )] } 


J= 1 

+ F 0 {r n <5 0n ±^< 1 »(2 t -A 0 /i)} 


n = 0, 1, • • ■ , oo, 


where, 


and, 


J n (zA 0 ) - A n J'(tA 0 ) 
(w - a ) 2 A 0 J„ (tAi) 


r„ = 


a„ = 


(3.5) 


(3.6) 


(3.7) 


Po^l^n^l) 

In the numerical calculations the series is truncated at s = N and then (3.5) yields 
a set of 7V + 1 homogeneous equations for F , . These may be written in matrix form, 


[A]F = 0, (3.8) 

where F is a vector of length TV + 1 of the unknown coefficients F , . For a non-trivial 
solution to exist the determinant of this matrix must be zero. This provides the 
dispersion relationship between the wavenumber and frequency. 
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A similar set of equations may be derived for the unknown coefficients in the 
series representing the solutions that are odd about the x — y plane. These may be 
written, 

F, = (-1)'F_,, (3.4) 


then the independent equations yielded by (3.2) may be written, 

E G - { r »' s - ± [«iV. (2*A„ A) + (-l)*l?y_ 1 .(2«-A 0 fe)] } 

9 = 1 

n = 0, 1, ■ • • , oo. 


(3.9) 


The requirement of a non-trivial solution for G t results in a dispersion relationship 
for the odd modes about the x — y plane. 


A similar expression to both eqns. (3.5) and (3.9) was obtained by SedeTnikov 
(1967). He developed dispersion relationships for multilayer jets, several jets, and 
jets between parallel walls or in rectangular ducts. In each case the jet was repre- 
sented by a vortex sheet. No roots of the dispersion relationship were determined. 
Written in the form of eqns. (3.5) and (3.9) the off-diagonal elements vanish for 
large jet separations and the eigenvalues are the zeroes of T n . These eigenvalues 
correspond to the axisymmetric and helical normal modes of a single jet. 


It is clear that this form of the equations does not hold for zero frequency. This 
case is of interest, as it may be used in a description of the shock cell structure of the 
jet. Tam and Tanna (1982) showed how a model for the shock cell structure could 
be posed as an initial-value problem in which the fully-expanded vortex sheet acts 
as a waveguide for the pressure perturbation at the jet exit. For the steady problem 
the matching conditions at the vortex sheet require that the pressure perturbation 
be zero at and outside the vortex sheet. Thus there is no communication between 
the two jets and the shock cell structure remains unchanged from the single jet case. 
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However, it should be noted that this is only true for jets into stationary air and 
some coupling between the steady shock cell structure could occur if the ambient 
air were in motion; see Morris(1988). 


3.2 Calculations 


There are many parameters and operating conditions that could be varied for 
the present configuration. Thus, the calculations have been limited to a set of 
operating conditions that correspond to available experiments. In the calculations 
for both the vortex sheet and the realistic mean profile representation of the jet 
the diameter of the jet is taken to be the fully- -expanded jet diameter. It may be 
argued, that in either the case of an over- or under-expanded jet, the jet plume 
will adjust its cross-section so as to preserve mass flux but equalize the mean static 
pressure. This gives the following relationship, assuming isentropic flow, between 
the fully-expanded and design jet dimensions. 




Mi 

Md 


1 + ^Ml 
1 + 


7 + l 
4 ( 7 - 1 ) 


(3.10) 


where r d is the non-dimensional design jet radius and M : and M d are the fully- 
expanded and design jet Mach numbers respectively. 


The instability wave calculations in the vortex sheet case provide an indication 
of the character of the results to be expected in the more realistic calculations that 
include the effects of finite mixing layer thickness. Four types of solution may be 
classified as shown in Table 1. In addition, each solution may be classified by the 
azimuthal mode number it approaches as the jets move further apart. 


Figure 2 shows the variation of the axial growth rate, —kj as a function of 
the separation distance between the two jets’ centerlines h. The instability wave 
frequency is 1.0 in each case. This corresponds to a Strouhal number of 0.318, 
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( 1 / 7 r) , based on the fully-expanded jet diameter and velocity. The design and fully- 
expanded jet Mach numbers are 1.0 and 1.32 respectively and the jet is unheated. 
These conditions correspond to the experiments of Seiner, Manning and Ponton 
(1988). In the numerical evaluation of the dispersion relationships obtained from 
(3.5) and (3.9) a value of N = 5 was used. Calculations were also performed with 
N = 9 with no significant change in the calculations, even for small values of jet 
separation. For large separations the solutions approach those for the single jet and 
the growth rate increases with azimuthal mode number. Calculations are presented 
for mode numbers 0, 1 and 2. For a given mode number the most unstable mode 
type is a function of jet spacing. For example, consider the mode number 1. For 
h > 1.6 the type III mode is the most unstable, though its value does not differ 
greatly from the value at large h. For h < 1.6 the type I mode dominates. For 
the range of mode numbers and types considered, the mode number 1, type I mode 
appears to be most affected by small separations showing a large increase in axial 
growth rate. This mode is associated with a motion that is even about both x — y 
and x — z planes and is dominated by a pair of helical motions of opposite sense 
about each jet. This is the type B mode described by Seiner, Manning and Ponton 
(1988) which they found to be the dominant mode in their twin plume resonance 
experiments. 

The prediction that the stability of the (1,1) mode is affected strongly by the jet 
separation and the interchange of dominance between modes of different types as the 
jet separation changes is encouraging, as it provides qualitative agreement with the 
observations of Seiner, Manning and Ponton (1988) and Wlezian (1987). However 
these results should be treated with some caution as they are based on the vortex 
sheet model for the jet. The dominant or preferred mode of a real jet is determined 
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by the total growth of a given frequency disturbance through the developing shear 
layer. In the next section more realistic profiles are used to describe the jet flow. 

4. REALISTIC JET FLOW MODEL 

4.1 Analysis 

In this case the mean velocity and density vary in Region II in a smooth, real- 
istic manner. The matching between the potential core and ambient flow solutions 
must be performed using a numerical solution in the mixing layer. 

In the potential core the solution for the pressure fluctuation takes the form 
given by (2.11). Consider, for example the modes that are even about the x — y 
plane. For n = — oo, • • • , oo (2.2) may be integrated from r = Ri to r = R 2 with 
initial conditions, 

P = J n (*Aj Ri ) and ~ = iX 1 J' n (iX 1 R 1 ). (4.1) 

The corresponding numerical solutions at r = R? are denoted by p n and p' n . These 
solutions may be matched with the exterior solutions for all n. That is, 

OO 

A n p n = ]T F,0. n ( 1 A 0 R 2 ), (4.2) 

8 = — OO 

and 

00 

1 = — OO 

Following the same approach used in Section 3, based on the symmetry properties 
of p n , /3 sn and F s , a dispersion relationship may be derived from an identical system 
of equations to (3.5). However, in this case the T n and A n are defined by, 

_ H^jiXoR*) - k n H l n 1) '{i\ 0 R*) (44) 

J n (fA 0 R 2 )-A n J'(t-A 0J R 2 ) 
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and, 


A„ = iX 0 p n /p' n . (4.5) 

For inodes that are odd about the x — y plane the dispersion relationship may be 
obtained from (3.9) with T„ and A n defined by eqns.(4.4) and (4.5) respectively. 


4.2 Calculations 

In the subsequent calculations it is assumed that the mean velocity and density 
of the jet flows take the same form as in the single jet case: up to the location where 
the jet edges meet. The mean velocity is assumed to take the form: 

X _ / 1 r<g{x) Z 4 6 \ 

u ( r > 1 ) | exp [- ln(2)rj 2 ] r > g{x). 


where, 

h = \r ~ g{x)]/b(x). (4-7) 


g(x) is the radius of the potential core and b( x) is the half-width of the mixing 
layer. 


The mean density is related to the mean velocity through a Crocco relationship, 


P = 


^ — ^u(l - u)M? + u + To (1 - «) 


(4.8) 


where T 0 is the non-dimensional ambient temperature. 

From the mean axial momentum integral equation a relationship may be found 
between the potential core radius and the half-width of the mixing layer, 

g{x) = -ft b + sJbHfil -2/? 2 ) + l, (4.9) 


where, 


Pi = f p u 2dr i, 

Jo 
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and 


/?a = / P u 2 r]drj. 

Jo 

At some axial location the edges of the two jets will touch on the symmetry 
plane and the present analysis, that assumes that the mean flow is axisymmetric 
relative to each jet’s centerline, is no longer valid. In the present calculations the 
edge of the jet is taken to be the location at which the axial velocity given by (4.6) 
equals 0.01. This corresponds to a value of tj of 2.58. Thus the pr< ent calculations 
are for values of jet thickness such that, 

g(b) + 2.586 < h. ( 4 - 10 ) 

A variable step-size fourth-order Runge-Kutta algorithm is used to integrate 
(2.2) from the edge of the potential core to r? = 2.58. This gives the values of p„ and 
p' n . As in the vortex sheet calculations the upper limit in the series representations 
is taken to be N = 5. Calculations have also been performed with N — 9 with 
negligible change in the largest elements of F, or G, . 

The vortex sheet calculations, shown in Fig. 2, indicate that for large separa- 
tions the higher azimuthal mode numbers have higher axial growth rates. Addi- 
tional calculations show that, for the present operating conditions, the maximum 
axial growth rate occurs for n — 3. A similar result is obtained for small values 
of local thickness 6(x). However, as the jet mixing layer thickens, the higher order 
azimuthal modes become damped more quickly. Figure 3 shows the variation of 
the axial growth rate —kj as a function of thickness 6(i) for a large jet separation 
h/r d — 5.0, for the first three azimuthal modes. In this and subsequent calculations 
the fully— expanded jet Mach number is 1.32, the design jet Mach number is 1.0, and 
the jet is unheated. Figure 3 shows how the growth rate of the n = 2 mode rapidly 
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decreases. The helical mode n = 1 has a larger growth rate than the axisymmetric 
mode n = 0 for the all values of jet thickness considered. In the subsequent calcula- 
tions only the two lowest mode numbers will be examined. The calculations shown 
in Figure 3 give results that are identical to the single jet case. 

As the separation decreases so the growth rates of the various mode numbers 
and types move away from their large separation value. Figure 4 shows this variation 
for mode numbers 0 and 1 and the four mode types. The jet thickness is b = 0.2 
and the Strouhal number St = 0.3. The change in the axial growth rate is relatively 
small for jet separations greater than 2 radii. For these conditions the most unstable 
mode at the closest separation achievable, before the jets’ edges merge, is the (1,IV) 
mode. This mode is dominated by two helical instabilities that are out of phase 
that give a solution that is odd about both the x — y and x — z planes. However, 
at other separations other modes are the most unstable. 

The relative instability of the various modes at a given separation has been 
found to be nearly independent of jet mixing layer thickness. For example, Figure 
5 shows the variation of axial growth rate with b(x) for the (0,1) and (0,111) modes. 
The single jet, n = 0 value is shown for comparison. In this case, with h/r d = 1.9, 
the (0,111) mode is the most unstable at all jet thicknesses. 

Before considering the eigenfunctions for the various modes of instability the 
effect of wave frequency will be considered. Figure 6 shows the variation of axial 
growth rate — /c ; with Strouhal number for the same modes shown in Fig. 5. The 
mixing layer thickness b = 0.2. Except for the lower Strouhal numbers the (0,111) 
mode is more unstable than the (0,1) mode or the n = 0 single jet mode. At 
this jet thickness the most unstable frequency occurs for a Strouhal number of 
approximately 0.45 and is relatively independent of jet separation or mode type. 
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The real part of the wavenumber is also affected by the jet separation. The 
trend in all the cases considered involves an increase in k R for the even modes 
about the x - z plane and a decrease in k R for the odd modes as the jet separation 
decreases. However, the changes are relatively small involving typically a 10% 
change from the single jet value. For example, Figure 7 shows the variation with 
mixing layer thickness of the phase velocity, given by cj/Zc^r, for the same modes 
shown in Figure 5. The phase velocity for the (0,1) mode, that is even about the 
x — z plane is lower than the single jet or large separation value. Conversely, the 
(0,111) mode, that is odd about the x - z plane takes a higher value. It should 
be noticed that for larger thicknesses, where the instability wave is reaching its 
maximum amplitude or neutrally stable condition, there is effectively no change in 
the phase velocity. In this region the phase velocity is approximately 0.73. Thus the 
observed shift in the screech frequency for the twin jets is linked to a change in the 
shock cell spacing rather than a modification to the phase velocity of the large scale 
structures. Seiner, Manning and Ponton (1988) did observe a 10 to 15% increase 
in the shock cell spacing. The reason for this increase is unclear as, in the absence 
of ambient flow, the shock cell structure of the two jets should be independent. 
However, insufficient aerodynamic data for the twin jets is available at present to 
help to explain this observation. 

The pressure distributions associated with each mode of instability may be 
constructed by obtaining the coefficients F t or G t for a given eigenvalue. An inverse 
iteration technique is used to obtain these values. That is, 

[A]F fc + 1 = <tF\ (4.11) 

where cr is a scaling factor and [A] is given by (3.8). An initial guess for F° is 
taken to be {1, 1, • ■ « , l} r . This algorithm has been found to give convergence in 


16 



two iterations for the cases considered. The interior coefficients may then be found 
from (4.2). 

For convenience in the present calculations only the pressure field outside the 
edge of the jet has been determined. Equation (2.3), rewritten in terms of modes 
that are odd or even about the x — y plane, is used to calculate the pressure. For 
example, Figure 8 shows contours of equal pressure level for the (0,1) mode and 
h/r d — 1.9. The phase, given by kx — u >t in (2.4) has been set to zero. It can be 
seen that the pressure field remains nearly axisymmetric. However, in the region 
between the two jets there is a loss of axisymmetry. In this region the amplitude 
of the pressure is nearly uniform and equal to the maximum amplitude achieved at 
the edge of the jet. The shaded region shows the region of maximum amplitude. 
This is the case for all the modes that are even about the x — z plane. 

A measure of the azimuthal mode content for each mode of instability and type 
is given by the relative magnitudes of the coefficients F, and G, . Table 2 shows how 
these amplitudes vary with jet separation for the ( 1 ,111) mode. For each separation 
the n — 1 helical mode is dominant. However, for hjr d — 1.5 the n = 0 mode 
amplitude rises to 67% of that of the n = 1 mode and the n = 2 mode rises to 16% 
of the n = 1 mode. 

5. DISCUSSION 

The present calculations have shown how the growth rates of instability waves 
or large structures in the initial mixing region of twin supersonic jets are affected by 
the jet separation. This interaction is caused by a coupling of the waves unsteady 
flow fields even before the time-averaged jet flows have merged. At a given operating 
condition the mode number and type that is most unstable is a function of the jet 
separation. However, the quantitative change in the local growth rates are relatively 
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small until the separation between the jet centerlines approaches the jet diameter: 
that is hjr d — * 1. Though this may be achieved in theory in the case of the vortex 
sheet model, similar calculations for the realistic jet profiles are limited to very 
small mixing layer thicknesses. At larger thicknesses the jets merge rendering the 
present analysis invalid. 

The influence of the jet separation on the eigenvalues can be seen from (3.5) to 
(3.8). The elements of the matrix A consist of two components: those that depend 
on /i, the jet separation, and those that do not. In fact, the easiest way to construct 
the matrix elements numerically is to first note that the components that depend 
on h form a symmetric matrix. The other terms, that only occur on the diagonal, 
may then be added. As the separation increases the relative magnitude of the terms 
that depend on h decrease rapidly, as Hankel functions, relative to the remaining 
components. Eventually, only the terms given by T n , that occur on the diagonal, are 
significant. The numerator in (3.6) can be seen to be the dispersion relationship for 
azimuthal mode number n for a single jet. Thus the single jet eigenvalues constitute 
the zeroes of the determinant of matrix A for large jet separations. 

It should be noted that the amplitude achieved by an instability wave depends 
on the integrated growth of the wave with axial distance and the local variation in 
the shape of the eigenfunction. Thus, relatively small changes in the local growth 
rate can result in large changes in the eventual amplitude of the wave. For example, 
using the data shown in Fig. 5, and assuming that dbjdx is given by the single jet 
value for the same operating conditions, the amplitude of the (0,111) mode is 24% 
larger for h/r d = 1.9 compared to h/r d — 5.0 at the location where the jets merge 
in the former case. However, as mentioned earlier, the rate of spread in the twin jet 
case may be decreased by the co-flowing, entrained air between the jets. This would 
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increase the relative amplitude at the merger location. In addition, the pressure 
levels between the jets are much higher in the twin jet case for modes I and II as 
shown in Fig. 8. However, it is not clear whether changes of this order of magnitude 
would be sufficient to explain the observed changes in the near-field pressure levels 
when the twin jets resonate. 

The present analysis has considered only a part of the feedback cycle associ- 
ated with twin jet screech. The instability wave’s growth into the merged jet region 
must be determined. In this region the merged jet would resemble more closely 
a developing rectangular jet. In this case certain normal modes, particularly the 
flapping mode about the x — y plane might be enhanced. The interaction between 
the instability waves and the shock cell structure in the jet that gives rise to the 
upstream propagating acoustic wave must then be described. It should be noted 
that existing theories of jet screech are unable to predict the amplitude even for 
single circular jets. Experimentally, the occurrence and amplitude of jet screech are 
very sensitive to small changes in the detailed geometry of the jet model and labo- 
ratory. So a prediction of the occurrence of resonance or its amplitude is extremely 
difficult. 

To assist in the extension of the present calculations to other sections of the 
feedback loop further experimental data on the aerodynamic development of the 
twin jets is required. This includes the modification to the rate of growth of the jet 
mixing layers, mean flow contours in the merged jet region, and measurements of 
the entrained flow between the jets. 

Though the present analysis does not answer all the questions regarding the 
complex phenomenon of twin jet resonance, it has shown how an instability wave 
analysis can provide some insight into the interaction of twin supersonic jets. 
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x — y plane 

x — z plane 

Mode 

even 

even 

I 

odd 

even 

II 

even 

odd 

III 

odd 

odd 

IV 

Table 1. Classification of Normal Modes. 


Mode 

number 

Relative amplitude 



hlr d = 3.0 

h/r d = 2.0 

h/r d = 1.5 

0 

0.016 

0.199 

0.199 

1 

1.000 

1.000 

1.000 

2 

0.001 

0.098 

0.155 

3 

0.000 

0.008 

0.019 


Table 2. Variation in relative mode amplitude with jet separation, Mode (l,III). 
M i - 1.32, M d = 1.0, St = 0.3, b = 0.2. 
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Figure Captions 


'Fig. 1 Schematic of twin jet cross section and coordinate systems. 

Fig. 2 Variation of the axial growth rate -kj with jet spacing. My = 1.32, M d = 

1.0, St = 0.318. , Mode I; , Mode II; , Mode III; Mode 

IV. 

Fig. 3 Variation of the axial growth rate — kj with jet mixing layer thickness. 

M y = 1.32, M d = 1.0, St = 0.3 ,h/r d = 5.0. , n=0; , n=l; 

> n=2. 

Fig. 4 Variation of the axial growth rate —kj with jet spacing. My = 1.32, M d — 

1.0, St = 0.3, 6 = 0.2. , Mode I; , Mode II; , Mode III; 

Mode IV. 

Fig. 5 Variation of the axial growth rate — kj with jet mixing layer thickness. 

My = 1.32, M d = 1.0, St = 0.3. , n = 0, h/r d = 5.0; , (0,1), 

h/r d — 1.9; , (0,111), h/r d = 1.9. 

Fig. 6 Variation of the axial growth rate —kj with Strouhal number. My = 
1.32, M d = 1.0, b = 0.2. For legend see Fig. 5. 

Fig. 7 Variation of the phase velocity, u/k R , with jet mixing layer thickness. My = 

1.32, M d = 1.0, St = 0.3. , n = 0, h/r d = 5.0; , (0,1), = 1.9; 

, (0,111), hjr d — 1.9. 

Fig. 8 Contours of equal pressure level, Mode (0,1). My = 1.32, M d = 1.0, St = 

0.3, h/r d = 1.9. , outer edge of the jet. Contours from -.025 to -.175 in steps 

of -.025. 
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Fig. 1 Schematic of twin jet cross section and coordinate systems. 
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